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SUMMARY 

PLRAY  (Ray  Progation  Loss)  is  a locally  generated  acoustic  ray  program 
similar  to  the  FACT  (Fast  Asymptotic  Coherent  Transmission)  model.  It  is  de- 
signed to  compute  propagation  loss  as  a function  of  receiver  range  at  constant 
depth  in  a horizontally  stratified  ocean  with  a flat  bottom.  Since  the  propa- 
gation loss  is  derived  from  the  resultant  intensity  of  the  multipath  propaga- 
tion, it  is  necessary  to  compute  the  intensities  of  the  individual  contributing 
rays  which  arrive  at  each  receiver  location.  As  in  FACT,  this  is  done  by  com- 
puting a family  of  rays,  specified  by  their  source  angles,  and  interpolating 
between  them  to  determine  the  intensities  of  those  rays  which,  if  they  had  been 
computed,  would  have  arrived  at  the  desired  receiver  location.  PLRAY,  like 
FACT,  offers  the  option  of  either  incoherent  or  semicoherent  ray  addition.  In 
the  incoherent  mode  of  summation,  the  individual  ray  intensities  are  added  di- 
rectly without  regard  to  phase.  In  the  semicoherent  mode  the  group  of  rays  in 
each  arrival  order,  such  as  the  four  single  bottom-bounce  rays,  are  treated 
coherently,  their  relative  phases  being  computed  in  a simple  approximate  manner. 
The  resultant  coherent  intensities  of  the  various  arrival  orders  are  then  added 
incoherently. 

PLRAY  is  similar  to  FACT  in  its  basic  structure.  It  processes  only  one 
source/receiver  depth  combination  at  a time,  but  computations  may  be  made  si- 
multaneously at  up  to  six  frequencies.  The  maximum  number  of  receiver  range 
points  is  251  in  PLRAY  and  250  in  FACT.  The  choice  between  incoherent  and 
semicoherent  ray  addition  may  be  made  separately  for  each  frequency. 

PLRAY  incorporates  wave  corrections  for  smooth  caustics  similarly  to  FACT 
but  lacks  a correction  for  cusped  caustics.  Although  the  absence  of  the  cusped 
caustic  correction  leads  to  excessively  large  errors  when  the  source  and  re- 
ceiver are  at  the  same  depth,  reasonable  results  are  obtained  when  the  source 
and  receiver  depths  differ  by  only  a few  feet.  When  the  source  and  receiver 
are  specified  at  or  very  near  the  same  depth,  the  program  automatically  shifts 
them  by  an  amount  necessary  to  maintain  a difference  of  two  feet. 

PLRAY  contains  a surface  duct  model,  similar  in  some  respects  to  that  of 
FACT,  which  is  activated  in  place  of  ray  computations  for  propagation  in  a sur- 
face duct.  The  surface  duct  model  is  used  when  both  the  source  and  receiver 
are  located  in  the  duct  or  when  one  is  in  the  duct  and  the  other  is  below  the 
duct  but  not  deeper  than  1.8  times  the  duct  depth. 

PLRAY  differs  from  FACT  in  the  manner  of  interpolating  between  the  points 
of  the  input  velocity  profile  table.  FACT  uses  linear  segments  to  connect  the 
points,  while  PLRAY  uses  curvilinear  segments  in  which  the  square  of  the  index 
of  refraction  (or  l/c2)  is  parabolic.  The  segments  are  joined  with  continuity 
of  slope.  PLRAY  also  differs  from  FACT  in  the  method  of  computing  ray  inten- 
sities. The  key  factor  in  the  intensity  computations  is  the  derivative  of  the 
range  with  respect  to  the  ray  source  angle  (or  some  monotonic  function  of  ray 
source  angle).  PLRAY  computes  the  derivative  from  analytical  formulas.  FACT 
fits  the  curve  of  range  vs  (a  function  of)  source  angle  with  a parabola  fitted 
to  three  points.  Such  a procedure  appears  to  be  advantageous  for  the  caustic 
correction  procedure  when  caustics  occur,  but  it  is  subject  to  inaccuracies  in 
the  computation  of  ordinary  ray  intensities. 
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PLRAY  also  differs  from  FACT  with  respect  to  bottom  loss  curves.  Both 
programs  contain  internally  the  sets  of  FNWC  (Fleet  Numerical  Weather  Central) 
curves.  FACT  contains  an  updated  set  of  nine  curves  for  the  frequency  3500  Hz, 
which  have  not  been  inserted  into  PLRAY.  The  programs  differ  in  the  manner  of 
implementing  the  curves  at  frequencies  other  than  the  specific  frequencies  for 
which  they  were  originally  devised.  PLRAY  interpolates  linearly  between  the 
pair  of  curves  on  either  side  of  the  frequency  of  interest.  FACT  does  not  in- 
terpolate but  divides  the  spectrum  into  frequency  intervals,  such  that  the 
bottom  loss  is  constant  throughout  each  interval  and  jumps  discontinuously  from 
one  interval  to  the  next.  With  PLRAY  the  user  is  not  limited  to  the  internally 
contained  FNWC  curves,  but  is  permitted  to  read  in  his  own  bottom  loss  curves 
in  their  stead. 

PLRAY  also  incorporates  a provision  for  reading  in  a beam  pattern.  The 
pattern  is  applied  at  the  source. 

PLRAY  has  been  compared  with  FACT  in  21  cases  involving  eight  different 
velocity  profiles.  In  most  of  the  cases  runs  were  made  both  semicoherently 
and  incoherently.  As  a check  on  the  accuracy  of  both  programs,  the  outputs  of 
the  semicoherent  runs  were  compared  with  the  output  of  a locally  generated  nor- 
mal mode  program  called  AP2.  To  facilitate  comparison  with  the  ray  program 
outputs,  the  normal  mode  curves  were  smoothed  by  application  of  a weighted 
sliding  window  which  removed  most  of  the  fine  structure. 

The  interference  patterns  generated  by  PLRAY  and  FACT  in  the  semicoherent 
runs  are  observed  chiefly  in  the  bottom-bounce  regions.  In  most  of  the  runs, 
computations  were  carried  out  sufficiently  far  in  range  to  cover  the  first  bot- 
tom-bounce region  and  a portion  of  the  second.  In  general,  both  ray  programs 
were  able  to  duplicate  in  remarkable  detail  the  features  of  the  AP2  normal  mode 
interference  pattern  in  the  first  bottom- bounce  region.  The  best  agreement  was 
observed  for  the  steeper  rays  at  short  ranges.  The  agreement  with  AP2  was  ob- 
served to  deteriorate  with  increasing  range,  increasing  receiver  (or  source) 
depth,  and  increasing  frequency.  With  a few  exceptions  the  agreement  was  rather 
poor  in  the  second  bottom-bounce  region.  Generally  poor  agreement  was  also  ob- 
served in  the  two  cases  in  which  both  source  and  receiver  were  located  within 
a depressed  sound  channel,  where  presumably  the  propagation  of  energy  in  the 
channel  tended  to  overshadow  the  bottom-bounce  propagation. 

In  many  of  the  runs  it  was  observed  that  while  the  three  programs  agreed 
quite  well  on  the  relative  locations  and  the  db  levels  of  the  individual  fea- 
tures of  the  interference  pattern,  there  was  a discrepancy  in  the  absolute 
ranges  appearing  as  an  expansion  or  contraction  of  the  range  scale.  The  dis- 
tortions are  usually  quite  small  at  short  ranges  but  tend  to  become  serious  at 
longer  ranges.  A number  of  possible  causes  of  the  scale  distortion  have  been 
examined,  but  no  completely  satisfactory  explanation  has  as  yet  been  found. 

Both  PLRAY  and  FACT  contain  algorithms  for  cutting  down  the  amplitude  of 
the  interference  oscillations  in  cases  where  the  number  of  receiver  range 
points  per  cycle  of  the  oscillations  is  insufficient  to  provide  adequate  sam- 
pling. The  coherence  factor  which  generates  the  oscillations  is  multiplied  by 
an  attenuation  coefficient  which  varies  with  the  number  of  points  per  cycle 
from  a full  value  of  1 when  the  sampling  is  adequate  to  0 when  the  sampling  is 
totally  inadequate.  PLRAY  and  FACT  differ  in  the  manner  in  which  the  number 
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of  points  per  cycle  is  computed.  PLRAY  computes  an  instantaneous  value  based 
on  an  analytical  formula,  while  FACT  computes  an  average  value  over  the  entire 
range  interval  covered  by  the  family  of  rays  involved.  Since  the  period  of  the 
oscillations  tends  to  vary  strongly  with  range  in  the  first  bottom-bounce  re- 
gion, cases  have  been  encountered  in  which  the  PLRAY  and  FACT  propagation  loss 
curves  differ  greatly  in  appearance;  with  the  PLRAY  curves  showing  a steady  in- 
crease in  amplitude  with  range,  and  the  FACT  curves  showing  a constant  ampli- 
tude, with  oscillations  exhibiting  a very  choppy  appearance  at  the  shorter 
ranges  due  to  inadequate  sampling.  However,  comparison  of  the  FACT  curves  with 
the  normal  mode  oscillations  indicates  that  adequate  sampling  can  be  obtained 
with  somewhat  fewer  points  per  cycle  than  are  called  for  in  the  algorithm. 
Further  investigation  of  the  sampling  problem  is  needed  with  a view  to  rede- 
fining the  parameters. 

Except  for  the  manner  in  which  the  sampling  problem  is  handled,  PLRAY  and 
FACT  yield  generally  comparable  results  in  the  semicoherent  runs  and  neither 
shows  any  clear  advantage  over  the  other. 

In  the  incoherent  runs  PLRAY  and  FACT  show  excellent  agreement  throughout 
the  bottom-bounce  regions  except  for  the  difference  in  the  ways  in  which  the 
internally  contained  bottom  loss  curves  are  applied. 

With  regard  to  the  prediction  of  convergence  zones,  PLRAY  and  FACT  agree 
with  AP2  on  the  locations  of  the  zones,  but  differ  widely  in  the  shapes  of  the 
zones  and  the  heights  of  the  peaks.  However,  even  though  the  caustic  correc- 
tions do  not  perform  as  well  as  might  be  desired,  they  are  a highly  significant 
improvement  over  unmodified  ray  theory. 

Runs  were  made  with  three  velocity  profiles  containing  surface  ducts  of 
appreciable  thickness.  In  those  runs  in  which  both  the  source  and  receiver 
were  located  within  the  duct  PLRAY  generally  agreed  closely  with  the  normal 
mode  predictions  in  the  range  interval  dominated  by  the  surface  duct  propaga- 
tion. The  results  of  FACT  were  somewhat  more  variable,  showing  good  agreement 
in  some  runs  and  less  satisfactory  agreement  in  others.  In  those  runs  in  which 
the  source  was  in  the  duct  and  the  receiver  slightly  below,  both  ray  programs 
gave  variable  results. 

In  the  runs  on  two  of  the  eight  velocity  profiles  FACT  produced  results 
of  questionable  validity.  One  of  the  two  profiles  contains  a shallow  SOFAR 
channel  with  a double  minimum.  The  source  and  receiver  were  located  on  oppo- 
site sides  of  the  first  minimum.  The  convergence  zones  predicted  by  FACT  have 
a strange  shape  consisting  on  the  near  side  of  a long  gradual  rise  covering 
about  20  kyd,  followed  on  the  far  side  by  a steep  drop.  The  other  profile 
contains  a large  surface  duct.  PLRAY  and  AP2  agreed  well  on  the  interference 
pattern  in  the  first  bottom-bounce  region.  FACT,  while  predicting  roughly  the 
correct  average  levels,  generated  an  entirely  difference  interference  pattern. 
When  the  receiver  was  placed  slightly  below  the  bottom  of  the  duct,  the  FACT 
program  failed  to  generate  any  interference  pattern  at  all  and  produced  virtu- 
ally identical  curves  at  frequencies  of  50  and  100  Hz.  Aside  from  the  cusped 
caustic  problem,  PLRAY  experienced  no  difficulties  of  comparable  magnitude  in 
any  of  the  runs. 
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CONCLUSIONS 

1.  PLRAY  and  FACT  in  general  yield  comparable  performance  in  runs  in 
deep  water  at  frequencies  below  1 kHz.  (No  comparisons  were  made  at  high  fre- 
quencies.) 


2.  The  interference  patterns  generated  by  both  PLRAY  and  FACT  in  the 
semicoherent  runs  show  good  agreement  in  the  first  bottom-bounce  region  with 
the  corresponding  patterns  generated  by  normal  mode  theory.  The  quality  of  the 
agreement  tends  to  deteriorate  with  increasing  range,  increasing  receiver  (or 
source)  depth,  and  increasing  frequency. 

3.  Distortions  of  the  interference  patterns,  which  appear  as  expansions 
or  contractions  of  the  range  scale,  are  frequently  observed  in  the  comparisons 
of  PLRAY,  FACT,  and  the  normal  mode  program.  These  distortions  do  not  normally 
affect  the  locations  of  convergence  zones. 

4.  Although  both  ray  programs  usually  agree  on  the  locations  of  conver- 
gence zones,  neither  PLRAY  nor  FACT  appears  capable  of  accurately  predicting 
the  shapes  of  the  zones  and  the  heights  of  the  zone  peaks  as  generated  by  nor- 
mal mode  theory.  However,  both  programs  offer  a great  improvement  over  conven- 
tional ray  theory,  which  predicts  highly  erroneous  spikes  at  caustics. 

5.  Except  for  differences  in  the  shapes  of  convergence  zones,  PLRAY  and 
FACT,  when  using  the  same  bottom  loss  curves,  produce  virtually  identical  re- 
sults in  the  incoherent  mode  of  ray  addition.  The  programs  differ,  however, 

in  the  manner  of  interpolating  with  respect  to  frequency  between  the  internally 
contained  FNWC  bottom  loss  curves.  (FACT  does  not  interpolate.) 

6.  The  parameters  in  the  algorithm  which  reduces  the  amplitude  of  the 
oscillations  of  the  interference  pattern  in  cases  of  inadequate  range  sampling 
are  not  properly  chosen.  Further  investigation  is  needed  with  a view  to  se- 
lecting optimum  values.  This  conclusion  applies  to  both  PLRAY  and  FACT. 

7.  The  surface  duct  model  of  PLRAY  yields  somewhat  more  accurate  results 
then  that  of  FACT  when  the  source  and  receiver  are  both  located  in  the  duct. 
Neither  program,  however,  proved  highly  accurate  when  the  source  was  in  the 
duct  and  the  receiver  was  slightly  below. 

8.  The  absence  of  a cusped  caustic  correction  in  PLRAY  leads  to  a gross 
error  when  the  source  and  receiver  are  at  the  same  depth.  However,  the  error 
is  reduced  to  an  acceptable  magnitude  with  negligible  effect  on  the  computa- 
tions at  other  ranges  simply  by  separating  the  source  and  receiver  by  a few 
feet  in  depth. 
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INTRODUCTION 

Approximately  10  years  ago  a computer  program  (reference  (a))  was  devel- 
oped at  the  Naval  Air  Development  Center  for  the  prediction  of  acoustic  propa- 
gation loss  in  a horizontally  stratified  ocean.  This  program  has  been  exten- 
sively used  over  the  past  decade  not  only  at  NAVAIRDEVCEN,  but  also  by  various 
other  naval  activities  and  private  contractors,  and  has  proved  to  be  a very 
valuable  tool.  However,  as  has  long  been  recognized,  it  suffers  from  a number 
of  limitations,  chief  among  which  are  the  following  : 

1.  Being  based  strictly  on  ray  theory,  it  is  subject  to  the  limitations 
of  ray  theory.  In  particular,  it  predicts  infinite  intensity  at  caustics  (ray 
focal  points),  and  it  fails  to  account  for  leakage  of  energy  by  diffraction 
out  of  sound  channels.  The  most  commonly  encountered  example  of  the  leakage 
problem  occurs  when  both  the  source  and  receiver  are  located  within  a surface 
duct.  In  such  cases  the  predictions  generated  by  the  program  may  be  grossly 
in  error,  particularly  at  low  frequencies. 

2.  Among  the  various  modes  of  operation  available  to  the  user  of  the  old 
NAVAIRDEVCEN  program,  the  so-called  CLC  option  was  designed  for  the  routine 
computation  of  propagation  loss  as  a function  of  range.  In  this  mode  of  oper- 
ation, the  program  computes  a large  family  of  rays  and  interpolates  between 
these  rays  to  compute  the  propagation  loss  at  the  specified  ranges.  The  CLC 
option  has  two  limitations.  First,  in  computing  the  multipath  propagation, 
only  incoherent  ray  summation  is  available;  that  is,  the  individual  ray  inten- 
sities are  added  directly  and  no  provision  is  made  for  including  the  effects 

of  phase  interference.  Secondly,  the  structure  of  the  program  and  the  size  of 
the  arrays  in  the  computer  are  such  that  rays  cannot  be  traced  beyond  the  range 
at  which  they  have  passed  15  vertices  (7-1/2  complete  ray  cycles).  Beyond  this 
point  the  rays  are  lost. 

3.  Coherent  ray  summation  is  available  when  the  target  ray  (eigenray) 
option  is  used.  In  this  mode  of  operation  each  specified  receiver  location  is 
treated  as  a separate  case,  and  a search  and  iteration  procedure  is  employed  to 
find  each  individual  ray  which  propagates  to  that  location.  The  principal 
limitation  of  the  target  ray  option  is  the  long  computer  running  time  required. 
An  additional  limitaticn,  which  has  proved  troublesome  in  a few  isolated  runs, 
is  the  restriction  of  the  number  of  target  rays  to  30. 

Although  the  second  and  third  of  the  above  limitations  could  be  removed 
by  re-structuring  the  program,  the  first  is  basic  because  it  is  inherent  in  the 
ray  approximation  to  wave  propagation.  Considerable  work  has  been  done  over 
the  past  decade  in  developing  wave  corrections  to  ray  theory  and  in  the  appli- 
cation of  tnese  corrections  to  practical  ray  programs.  The  best  known  and  most 
widely  used  of  these  programs  is  the  FACT  model  (reference  (b) ) , developed  by 
the  former  AESD  (Acoustic  Environmental  Support  Detachment),  Office  of  Naval 
Research.  This  model  incorporates  a number  of  new  features,  the  most  important 
of  which  are  caustic  corrections  and  an  approximate  method  of  computing  phase 
interference  effects  near  the  ocean  surface.  The  FACT  program  also  includes  a 
surface  duct  model  which  was  originally  developed  for  the  FNWC,  Monterey,  Cal- 
ifornia (reference  (c)),  and  a wave  correction  for  rays  propagating  in  a de- 
pressed sound  channel  when  the  source  and  receiver  are  both  located  near  the 
channel  axis. 
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The  FACT  program  has  been  made  available  to  the  Naval  Air  Development 
Center  and  is  currently  operational  on  the  NAVAIRDEVCEN  CDC  6600  computer. 
However,  after  an  examination  of  FACT  and  a brief  period  of  testing  it,  a de- 
cision was  made  to  develop  a somewhat  similar  program  locally.  One  considera- 
tion in  this  decision  was  a desire  to  retain  some  of  the  more  desirable  fea- 
ture of  the  old  NAVAIRDEVCEN  program,  especially  the  fitting  of  the  input  ve- 
locity profile  points  with  curvilinear  segments  and  the  attendant  procedure  for 
calculating  ray  paths  and  intensities.  The  FACT  program,  even  though  it  rep- 
resents a significant  advance  over  conventional  ray  programs  with  respect  to 
its  excellent  novel  features,  leaves  much  to  be  desired  in  those  portions  where 
standard  ray  theory  is  applied.  Instead  of  fitting  the  velocity  profile  with 
curvilinear  segments,  which  provide  continuity  of  slope,  it  employs  the  old 
method  of  straight  lines;  then  in  the  computation  of  ray  intensities  it  uses 
an  approximate  method  of  questionable  accuracy.  The  key  parameter  entering 
the  computation  is  the  derivative  of  the  ray  range  with  respect  to  the  ray 
source  angle  (or  some  monotonic  function  of  the  source  angle) . Instead  of 
computing  the  derivative  from  an  analytical  formula,  the  FACT  program  fits  a 
parabola  to  three  points  and  calculates  the  derivative  from  the  slope  of  the 
parabola.  Cases  have  been  found  where  this  method  leads  to  significant  errors. 

Another  consideration  in  the  development  of  a local  program  is  the  ease 
with  which  the  program  may  be  modified  to  meet  special  requirements  which  con- 
tinually arise  in  the  laboratory's  work.  It  is  far  simpler  (and  safer)  to  make 
internal  changes  in  a home-grown  program  whose  details  are  thoroughly  under- 
stood than  to  dig  through  the  intricacies  of  an  externally  prepared  program 
such  as  FACT. 

The  resulting  program,  called  PLRAY,  is  similar  to  FACT  in  its  basic 
structure  and  requires  essentially  the  same  computer  running  time.  In  any  one 
run  it  computes  the  propagation  loss  as  a function  of  range  for  a single  source/ 
receiver  depth  combination,  while  tip  to  six  frequencies  may  be  processed  simul- 
taneously. It  has  a basically  similar  coherence  feature,  so  that  the  propaga- 
tion loss  may  be  computed  either  incoherently  or  semicoherently . (Both  pro- 
grams have  a third  option  labeled  "coherent,"  but  this  is  not  truly  coherent  in 
either  program  and  appears  to  be  of  little  practical  value.)  PLRAY  incorpo- 
rates corrections  for  smooth  caustics  but  lacks  the  cusped  caustic  correction 
feature  of  FACT.  The  presence  of  cusped  caustics  leads  to  large  errors  when 
the  source  and  receiver  are  at  the  same  depth.  In  PLRAY  the  errors  are  reduced 
to  acceptable  proportions  by  moving  the  source  and  receiver  depths  slightly 
apart.  PLRAY  also  lacks  the  wave  correction  for  depressed  sound  channel  propa- 
gation, but  it  contains  a locally  generated  surface  duct  model  broadly  similar 
to  the  surface  duct  model  of  FACT.  Finally,  PLRAY  contains  a provision  for 
inserting  a set  of  beam  patterns. 
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GENERAL  DESCRIPTION  OF  THE  PROGRAM 

PLRAY  computes  propagation  loss  as  a function  of  range  at  constant  depth 
in  an  ideal  horizontally  stratified  ocean  whose  acoustic  structure  is  specified 
by  a single  velocity  profile.  It  does  this  by  computing  a family  of  rays  and 
interpolating  between  these  rays  to  obtain  the  intensities  at  the  various  spec- 
ified receiver  locations.  The  intensity  at  any  one  receiver  location  is  the 
result  of  the  contributions  of  a number  of  rays  which  propagate  over  different 
paths.  The  program  offers  the  option  of  selecting  either  of  two  different 
methods  of  computing  the  resultant  multipath  intensity:  (1)  incoherent  summa- 
tion, in  which  the  individual  ray  intensities  are  added  without  regard  to  phase 
interference,  and  (2)  semicoherent  summation,  in  which  the  group  of  rays  asso- 
ciated with  each  arrival  order  (e.g.,  the  four  single  bottom-bounce  rays)  are 
added  coherently  and  the  resultant  intensities  of  the  various  arrival  orders 
are  added  incoherently. 

The  program  is  written  in  FORTRAN  V3.0-P380  for  the  CDC  6600  computer. 
INPUTS 

The  first  data  card  contains  a set  of  integers,  including  the  number  of 
frequencies,  the  number  of  receiver  depths,  and  several  other  parameters  which 
are  discussed  in  the  detailed  description  of  the  program.  This  is  followed  by 
a set  of  environmental  input  cards.  These  inputs  consist  of  the  sea  state, 
bottom  parameter,  velocity  profile,  and,  if  desired,  a set  of  bottom  loss 
curves  which  the  user  may  wish  to  substitute  for  the  internally  contained 
curves.  The  environmental  inputs  are  followed  by  a set  of  up  to  six  frequen- 
cies, a single  source  depth,  a set  of  up  to  eight  receiver  depths,  and  a range 
card  which  contains  the  minimum  range,  range  increment,  and  maximum  range.  The 
minimum  range  (not  included  in  FACT)  permits  the  user  to  begin  the  computations 
at  any  range  desired.  The  maximum  number  of  range  points  allowed  is  251.  Fi- 
nally, if  desired,  a set  of  beam  patterns  may  be  inserted  as  tables  of  pressure 
ratio  (relative  to  the  beam  axis)  vs  ray  source  angle,  one  table  for  each  of 
the  frequencies  specified.  Beam  inputs  are  read  and  pattern  computations  are 
performed  in  SUBROUTINE  BEAM.* 

CURVE-FITTING  (SUBROUTINES  CURFIT,  PARAB,  FIT,  BRIDGE) 

Before  curve- fitting,  corrections  for  the  curvature  of  the  earth  are  ap- 
plied to  the  depths  and  sound  speeds  of  the  velocity  profile.  The  curve-fit- 
ting routine,  which  fits  curvilinear  segments  which  join  with  continuous  slope, 
is  taken  directly  from  reference  (a)  and  will  not  be  described  further  in  this 
report . 


* Experience  has  shown  that  the  above  method  of  specifying  the  beam  patterns 
is  cumbersome.  It  is  included  in  the  current  description  because  it  is  per- 
fectly general  and  does  not  contain  any  classified  information.  For  any 
given  application  it  is  usually  preferable  to  specify  the  physical  parame- 
ters of  the  beam  hardware  and  let  the  computer  program  calculate  the  beam 
response.  This  is  readily  accomplished  simply  by  rewriting  the  BEAM  sub- 
routine without  affecting  the  rest  of  the  program. 
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RECEIVER  DEPTH  LOOP 

At  this  point  a loop  over  receiver  depths  is  entered,  which  extends  to  the 
end  of  the  program.  All  operations  from  this  point  on  must  be  carried  out 
separately  for  each  different  receiver  depth.  At  the  beginning  of  the  loop 
the  source  and  receiver  depths  are  compared  to  determine  whether  they  differ  by 
at  least  2 feet.  If  not,  they  are  automatically  shifted  to  maintain  this  sep- 
aration. 

AUGMENTATION  OF  PROFILE  (SUBROUTINE  INSERT) 

Before  proceeding  with  the  ray  computations,  a new  velocity  profile  table 
is  set  up  by  inserting  the  source  and  receiver  depths  into  the  original  pro- 
file. The  sound  speeds  at  the  source  and  receiver  depths  are  calculated  from 
the  parameters  of  the  respective  curvilinear  segments  in  which  they  lie.  At 
this  point  the  source  and  receiver  are  redefined  for  the  purpose  of  the  sub- 
sequent computations.  The  sound  speeds  at  source  and  receiver  are  compared, 
and  the  depth  with  the  larger  sound  speed  is  treated  as  though  it  were  the 
source,  while  the  depth  with  the  smaller  sound  speed  is  treated  as  though  it 
were  the  receiver.  In  what  follows,  the  original  pair  of  points  will  be  la- 
beled the  true  source  and  true  receiver,  and  the  new  pair  will  be  labeled  the 
ray  source  and  ray  receiver,  or  simply  the  source  and  receiver. 

SECTORS  (SUBROUTINE  SECLIM) 

The  family  of  rays  computed  actually  consists  of  a group  of  sub-families 
of  rays  of  different  types.  The  sub-families  are  bounded  by  limiting  rays. 

For  example,  the  limiting  ray  to  the  ocean  bottom  separates  the  RSR  rays  from 
the  bottom-bounce  rays.  In  order  that  the  interpolations  may  be  performed,  it 
is  necessary  to  separate  the  sub-families  and  work  with  each  of  them  one  at  a 
time. 


For  this  purpose  it  is  convenient  to  specify  the  rays  in  terms  of  the  ab- 
solute values  of  their  source  angles.  Each  sub-family  is  then  confined  to  a 
sector  in  the  angular  space  from  0 to  90  degrees,  the  sectors  being  bounded  by 
the  appropriate  limiting  ray  source  angles.  The  sectors  are  processed  one  at 
a time  in  order  of  increasing  source  angles.  The  outermost  sector  is  terminated 
at  30  degrees.  No  rays  are  computed  at  steeper  angles.  Instead,  the  intensi- 
ties are  computed  by  a simplified  approximate  formula. 

RAY  SUB-FAMILIES  (SUBROUTINE  DIVSEC) 

The  number  and  spacing  of  the  rays  of  each  sub-family  are  automatically 
determined  by  an  algorithm  in  the  program.  The  spacing  is  extremely  close  near 
the  inner  edge  of  the  sector  and  increases  in  a geometric  progression  until  it 
reaches  a constant  limiting  value.  Normally  the  constant  spacing  is  continued 
until  the  outer  edge  is  reached.  However,  if  the  sector  is  bounded  at  its 
outer  edge  by  a limiting  ray  to  a local  maximum  of  the  velocity  profile  (as 
distinguished  from  a limiting  ray  to  the  surface  or  bottom) , the  spacing  is  de- 
creased, on  approaching  the  boundary,  in  a manner  symmetric  to  the  spacing  at 
the  inner  edge. 
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RAY  COMPUTATIONS  (SUBROUTINE  INCR) 

Three  basic  horizontal  range  intervals  are  computed  for  each  ray  of  the 
sub- family:  the  range  interval  from  source  depth  to  receiver  depth,  the  range 
interval  corresponding  to  one  complete  ray  cycle,  and  the  range  interval  cor- 
responding to  the  upper  portion  of  a cycle  above  the  receiver  depth.  These 
three  quantities  supply  all  the  information  necessary  to  compute  the  ray  ranges. 
At  the  same  time  the  corresponding  three  intervals  are  computed  for  the  range 
derivatives,  which  are  required  for  intensity  computations.  In  PLRAY,  the  de- 
rivatives are  taken  with  respect  to  the  tangent  of  the  source  angle. 

ARRIVAL  ORDERS 

Although  applicable  to  sources  and  receivers  at  any  depth,  the  structure 
of  both  FACT  and  PLRAY  with  respect  to  arrival  orders  is  most  easily  explained 
when  the  source  and  receiver  are  located  in  the  upper  portion  of  the  ocean.  In 
this  case  the  bulk  of  the  ray  cycle  consists  of  the  portion  which  passes  through 
the  deep  ocean,  and  the  direct  and  surface-reflected  (or  upper-refracted)  rays 
are  relatively  close  together. 

The  arrival  order  number  is  defined  as  the  number  of  passages  through  the 
deep  ocean,  that  is,  the  number  of  bottom  bounces  or  deep  refractions.  The 
zero  order  arrival  consists  of  the  two  rays  which  propagate  directly  to  the  re- 
ceiver without  going  deep.  The  first  and  subsequent  arrival  orders  each  con- 
sist of  four  rays.  If  the  source  and  receiver  are  shallow,  all  the  rays  in  a 
given  arrival  order  propagate  over  similar  paths  and  may  be  expected  to  exhibit 
some  degree  of  coherence.  If  the  source  is  deep  and  the  receiver  shallow,  co- 
herence is  to  be  expected  only  at  the  receiver  end,  in  which  case  the  four  rays 
are  treated  as  two  pairs,  one  pair  leaving  the  source  in  a downward  direction 
and  the  other  pair  in  an  upward  direction.  Reciprocal  considerations  apply  to 
the  case  where  the  source  is  shallow  and  the  receiver  is  deep.  Within  each 
sector  the  arrival  orders  are  processed  sequentially,  beginning  with  order  zero. 

The  individual  ray  types  (such  as  down-up,  up-down-up,  down-up-down,  and 
up-down-up-down)  which  comprise  the  multipath  system  of  any  given  arrival  order 
will  be  termed  the  multipath  ray  types. 

PROCESSING  OF  MULTIPLE  FREQUENCIES 

All  operations  up  to  this  point  have  been  independent  of  frequency  and  the 
results  thus  far  apply  equally  to  all  frequencies  which  have  been  specified. 
However,  from  this  point  on  the  factors  which  enter  the  computation  of  intensi- 
ties are  frequency  dependent.  In  order  that  all  the  frequencies  of  interest 
may  be  processed  simultaneously,  arrays  are  provided  as  necessary  in  the  pro- 
gram for  the  frequency-dependent  variables  and  frequency  DO  loops  are  provided 
as  necessary  for  the  computations.  It  is  to  be  understood  that  the  description 
which  follows  applies  to  each  of  the  frequencies  which  are  being  simultaneously 
processed . 

CAUSTICS  (SUBROUTINE  RCAUST) 

Before  ray  intensities  are  computed  for  any  given  arrival  order  within  any 
sector,  a search  is  conducted  (in  SUBROUTINE  RAYSUM)  to  determine  whether  any 
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caustics  exist.  The  range  vs  source  angle  curve  for  each  multipath  type  is  ex- 
amined to  see  whether  a minimum  or  maximum  occurs  within  the  sector.  Since  by 
ray  theory  the  intensity  is  inversely  proportional  to  the  range  derivative,  the 
predicted  intensity  becomes  infinite  at  a minimum  or  maximum  of  the  range  curve. 
This  is  the  location  of  a caustic. 

If  a caustic  is  found,  caustic-correction  intensities  based  on  the  Airy 
function  Ai  are  computed  over  a range  interval  extending  from  a predetermined 
limit  on  the  shadow-zone  side  of  the  caustic  to  a predetermined  limit  on  the 
illuminated  side.  If  a beam  has  been  specified,  the  intensities  at  all  re- 
ceiver range  points  in  this  range  interval  are  multiplied  by  the  square  of  the 
beam  pressure  ratio,  evaluated  at  the  source  angle  of  the  ray  to  the  caustic. 

The  range  interval  is  then  marked  with  a flag  so  that  it  may  be  excluded  from 
the  subsequent  ray  intensity  computations. 

RAY  INTENSITY  COMPUTATIONS  (SUBROUTINE  RAYSUM) 

After  the  caustic-correction  intensities  have  been  computed  in  the  pre- 
scribed range  interval,  ray  intensities  are  computed  at  the  remaining  receiver 
locations  covered  by  the  rays  of  the  sector.  If  a beam  has  been  specified,  the 
ray  intensities  are  multiplied  by  the  squares  of  the  respective  beam  pressure 
ratios.  At  each  receiver  location  the  sum  of  the  intensities  of  the  contrib- 
uting multipath  ray  types  is  calculated,  yielding  the  resultant  incoherent  in- 
tensity for  the  current  arrival  order.  If  incoherent  ray  summation  has  been 
called  for,  this  resultant  value  is  added  to  whatever  value  is  already  present 
in  the  pertinent  range -frequency  bin  of  the  resultant  intensity  array. 

COHERENCE  COMPUTATIONS  (SUBROUTINE  RAYSUM) 

If  semi-coherent  addition  is  called  for,  it  is  necessary  to  calculate  a 
coherence  factor  by  which  the  incoherent  intensity  sum  is  multiplied  before 
being  stored  in  the  resultant  intensity  array.  (There  are  special  cases  in 
which  the  above  procedure  is  slightly  modified;  see  detailed  description  of  the 
program.)  It  should  be  mentioned  at  the  outset  that  coherence  is  calculated 
only  for  those  rays  which  strike  the  surface. 

The  coherence  factor  is  computed  in  a simple  approximate  manner.  Consider 
a pair  of  rays  which  leave  the  source  at  approximately  the  same  angle,  one  go- 
ing upward  to  be  reflected  at  the  surface  and  the  other  going  directly  down- 
ward. After  the  surface  reflection,  both  of  these  rays  travel  along  essentially 
the  same  path  and  arrive  together  at  the  receiver.  The  reflected  ray  may  be 
considered  to  have  originated  at  an  image  source  at  a point  above  the  surface, 
located  symmetrically  with  respect  to  the  actual  source.  The  assumption  is 
made  that  in  the  vicinity  of  the  source  these  two  rays  are  parallel  straight 
lines  and  that  their  difference  in  phase  is  simply  the  result  of  the  difference 
in  path  length,  plus  a 180-degree  phase  shift  resulting  from  the  surface  re- 
flection. A similar  effect  occurs  at  the  receiver  end.  These  effects  are  for- 
mulated mathematically  in  a coherence  factor  which  multiplies  the  incoherent 
intensity  sum.  As  the  range  to  the  receiver  varies,  the  phase  differences  also 
vary,  leading  to  an  oscillatory  interference  pattern. 

Two  limitations  are  placed  on  the  coherence  factor.  First,  if  the  oscil- 
lations of  the  interference  pattern  occur  rapidly  in  range,  there  may  be  a 
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sampling  problem.  The  parameter  of  interest  here  is  the  number  of  specified 
receiver  range  points  per  cycle  of  oscillation  of  the  pattern.  If  the  number 
of  points  per  cycle  is  too  small,  then  when  the  propagation  loss  is  plotted 
against  range  the  resulting  curve  will  not  adequately  represent  the  situation. 
This  problem  is  treated  in  a somewhat  arbitrary  way  by  simply  reducing  the  am- 
plitude of  the  oscillations.  An  attenuation  factor  is  introduced,  which  is  a 
function  of  the  number  of  points  per  cycle. 

The  second  limitation  is  concerned  with  the  tendency  for  the  coherence  to 
deteriorate  when  the  two  interfering  rays  are  separated  too  far  in  range.  This 
situation  may  be  expected  to  occur  when  the  source  or  receiver  is  deep  and  will 
tend  to  be  aggravated  at  shallow  angles,  where  the  assumption  of  straight-line 
propagation  becomes  less  valid.  The  problem  is  handled  by  arbitrarily  assuming 
that  coherence  occurs  only  when  the  horizontal  separation  between  the  two  in- 
terfering rays  at  the  source  depth  (for  interference  at  the  source  end)  or  at 
the  receiver  depth  (for  interference  at  the  receiver  end)  is  less  than  a pre- 
determined limiting  value. 

SURFACE  DUCT  PROPAGATION  (SUBROUTINE  SDUCT) 

If  the  source  and  receiver  are  located  within  a surface  duct,  the  use  of 
ray  theory  to  calculate  the  propagation  loss  does  not  in  general  lead  to  satis- 
factory results.  The  greatest  source  of  error  in  shallow  ducts  and/or  at  low 
frequencies  is  the  leakage  of  energy  out  of  the  duct.  Leakage  is  not  accounted 
for  in  ray  theory.  In  addition,  even  in  good  ducts,  early  experience  with 
PLRAY  has  indicated  that  the  caustic  correction  procedure  cannot  adequately 
handle  the  many  near-horizontal  caustics  occurring  in  the  ducts.  For  these 
reasons  PLRAY  incorporates  a surface  duct  model  which  calculates  the  propaga- 
tion loss  designed  to  approximate  the  predictions  of  normal  mode  theory  for 
propagation  within  and  immediately  below  the  duct. 

PROPAGATION  LOSS  (PROGRAM  PLRAY) 

After  all  intensities  have  been  calculated  they  are  converted  to  propaga- 
tion loss  and  stored  again  in  the  same  array. 

PLOT  ROUTINE  (SUBROUTINE  PLPLOT) 

The  program  contains  a plot  routine  for  the  purpose  of  plotting  the  prop- 
agation loss  vs  range  on  a Houston  COMPLOT,  Model  DP-1. 

CORE  REQUIREMENT 

Without  the  software  for  the  Houston  plotter  PLRAY  requires  slightly  less 
than  21000  (decimal)  words  on  the  CDC  6600  computer.  The  Houston  plotter  adds 
approximately  1600  words. 

RUNNING  TIME 

PLRAY  and  FACT  require  essentially  the  same  running  time.  For  example, 
in  one  of  the  typical  runs  made  for  this  report  the  time  was  12.46  sec  for 
PLRAY  and  11.66  sec  for  FACT. 
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DETAILED  DESCRIPTION  OF  THE  PROGRAM 

As  an  aid  in  understanding  the  following  detailed  description,  the  reader 
is  referred  to  the  FORTRAN  listing  of  the  program  in  appendix  A and  to  the 
glossary  of  the  principal  FORTRAN  symbols  presented  in  appendix  B.  A thumb- 
nail description  of  the  various  subroutines  and  functions  is  given  in  appendix 
C.  A detailed  description  of  the  input  data  deck  for  a run  of  PLRAY  will  be 
found  in  appendix  D. 

DATA  INPUT  SECTION 

The  first  portion  of  PLRAY  is  devoted  principally  to  reading  and  printing 
the  input  data,  applying  corrections  for  the  earth's  curvature  to  depths  and 
sound  speeds,  and  fitting  the  velocity  profile  with  curvilinear  segments. 

Card  1 contains  7 integers.  NF  is  the  number  of  frequencies  for  which 
computations  are  to  be  run,  not  to  exceed  6.  It  also  performs  the  additional 
function  of  serving  as  a code  to  stop  the  run.  The  program  is  written  in  the 
form  of  a huge  loop.  After  the  computations  have  been  completed  for  a set  of 
input  data,  the  program  returns  to  the  first  READ  statement  and  checks  the 
variable  NF.  If  NF  has  a nonzero  value,  a new  set  of  input  data  is  read  in  and 
processed.  However,  if  NF  is  zero,  the  program  stops;  thus,  a blank  card  may 
be  used  to  stop. 

NREC  is  the  number  of  receiver  depths,  limited  to  8. 

IPROF  is  a code  parameter  governing  the  environmental  inputs,  namely,  the 
velocity  profile  and  the  externally  specified  bottom  loss  curves.  A value  of 
zero  is  selected  for  normal  operation.  If  IPROF  is  set  to  1,  the  READ  state- 
ments for  the  velocity  profile  table  and  the  bottom  loss  tables  are  bypassed. 

This  value  is  selected  when  it  is  desired  to  run  an  additional  set  of  computa- 

tions using  the  same  environmental  inputs  (for  example,  an  additional  source 
depth  or  frequency  or  a different  set  of  ranges).  When  IPROF  is  set  to  2,  only 
the  velocity  profile  inputs  are  bypassed,  and  when  set  to  3,  only  the  bottom 
loss  inputs  are  bypassed.  (It  should  be  pointed  out  that  the  bottom  loss  in- 
puts are  bypassed  also  for  all  values  of  the  bottom  parameter  IBP  except  6.) 

The  integer  IUNIT  is  required  when  the  velocity  profile  is  specified  in 
terms  of  temperatures  and  salinities  instead  of  sound  speeds.  If  IUNIT  = 0, 
it  is  assumed  that  the  temperatures  are  in  °F;  if  IUNIT  = 1,  they  are  in  °C. 

I PLOT  is  the  plot  control  parameter.  If  IPLOT  is  assigned  a value  of  1, 

a plot  is  called  for. 

IPR  is  a print  control  parameter.  For  normal  output,  IPR  is  set  to  zero 
(or  left  blank).  A value  of  1 causes  the  program  to  print  a certain  amount  of 
diagnostic  information  which  is  useful  in  analyzing  the  details  of  the  run.  A 
value  of  2 produces  more  detailed  diagnostic  information. 

IBEAM  is  the  beam  control  parameter.  The  default  value  0 signifies  that 
there  is  no  beam  pattern.  If  IBEAM  = 1,  a set  of  NF  beam  patterns  is  to  be 
read  in.  If  an  additional  run  is  to  be  made  using  the  same  beam  patterns, 
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IBEAM  is  set  to  2.  In  this  case  the  beam  inputs  are  bypassed  and  the  previ- 
ously inserted  beam  patterns  apply. 

Card  2 is  the  environmental  input  header  card.  The  first  variable,  ISP, 
is  the  sea  state.  It  is  used  in  computing  the  surface  reflection  loss.  IBP 
is  the  bottom  parameter.  Values  from  1 through  5 call  for  application  of  the 
internally  contained  FNWC  bottom  loss  curves,  types  1 through  5.  The  exter- 
nal bottom  loss  input  statements  are  bypassed.  A value  of  6 signifies  that  ex- 
ternal bottom  loss  data  are  to  be  read  in.  LAT  is  the  latitude  of  the  velocity 
profile,  a real  (decimal)  number.  The  latitude  is  required  for  the  computation 
of  sound  speed  from  temperature  and  salinity.  The  default  value  is  45  degrees. 
When  the  sound  speeds  are  read  in  directly,  LAT  serves  only  for  identification. 
IDR  is  a letter  N or  S and  serves  only  for  identification.  IDBT  is  a run  iden- 
tification parameter  consisting  of  up  to  40  alphanumeric  characters. 

Card  set  3 is  a set  of  cards  on  which  the  external  bottom  loss  curves  are 
specified.  These  cards  are  needed  only  when  IBP  = 6 and  IPROF  = 0 or  2.  The 
first  of  these,  card  3a,  contains  the  variables  NFR,  the  number  of  frequencies 
for  which  curves  are  to  be  provided,  and  FR(J),  the  values  of  these  frequencies 
in  Hz.  The  FR(J)  must  be  arranged  in  a monotonically  increasing  order.  The 
bottom  loss  values  in  db  for  each  frequency  are  specified  on  a pair  of  cards 
3b.  Fifteen  points  are  required.  The  last  value  specified  must  correspond  to 
90  degrees.  If  less  than  15  points  are  required  to  specify  the  curve,  the  re- 
maining 10-column  fields  may  be  left  blank.  The  corresponding  grazing  angles 
in  degrees  are  specified  on  a pair  of  cards  3c.  The  cards  are  stacked  in  the 
order  3b,  3c  for  the  first  frequency,  3b,  3c  for  the  second  frequency,  etc. 

Card  set  4 is  the  set  of  cards  containing  the  velocity  profile,  one  caTd 
for  each  point.  Each  card  contains  fields  for  the  depth  ZB,  sound  speed  CB, 
temperature  T,  and  salinity  S.  If  sound  speeds  are  to  be  read  in,  the  temper- 
ature and  salinity  are  not  required.  If  temperature  and  salinity  are  to  be 
read  in,  the  sound  speed  field  must  be  left  blank  (or  zero).  The  integer  NDCD 
is  a code  variable  which  terminates  the  read-in.  A nonzero  value  must  be  in- 
serted for  NDCD  on  the  last  card.  These  cards  must  be  stacked  in  order  of  in- 
creasing depth. 

The  profile  variables  may  be  specified  in  either  English  units  (ft,  ft/sec, 
°F)  or  metric  units  (m,  m/sec,  °C) . If  sound  speeds  are  read  in,  the  program 
automatically  distinguishes  between  the  two  systems  of  units.  If  temperatures 
are  read  in,  IUNIT  must  be  0 for  the  English  system  and  1 for  the  metric  sys- 
tem. When  the  inputs  are  temperatures  and  salinities,  the  sound  speeds  are 
calculated  from  Wilson's  equation  (reference  (d))  in  FUNCTION  VELOC. 

As  each  depth  and  sound  speed  are  read  in  (or  calculated) , corrections  are 
applied  for  the  curvature  of  the  earth  (see  appendix  E) . Upon  completion  of 
the  profile  input,  SUBROUTINE  CURFIT  is  called  for  the  purpose  of  fitting  the 
profile  with  curvilinear  segments.  CURFIT  and  its  associated  subroutines  PARAB, 
FIT,  and  BRIDGE  are  taken  directly  from  the  old  NAVAIRDEVCEN  program  and  are 
described  in  reference  (a) . 

After  the  curve  fit  a test  is  made  for  the  presence  of  a surface  duct.  A 
search  is  made  for  the  largest  sound  speed  in  the  first  1000  feet  below  the 
surface  of  the  ocean,  and  a test  is  made  to  determine  if  this  largest  value 
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occurs  at  a local  maximum  point  of  the  profile.  If  so,  the  value  of  the  index 
IB  at  this  depth  is  stored  as  the  duct  parameter  IDUCT.  The  purpose  of  the 
test  is  to  determine  the  duct  depth  for  use  in  the  surface  duct  model.  If  the 
maximum  sound  speed  occurs  at  the  surface,  there  is  of  course  no  surface  duct, 
and  IDUCT  =1.  If  a duct  is  observed  but  is  found  to  be  deeper  than  1000  ft, 
IDUCT  is  set  to  1 and  the  surface  duct  model  is  not  implemented. 

Card  5 contains  the  NF  frequencies  F(IF)  and  the  corresponding  NF  values 
of  the  coherence  parameter  JCOH(IF) , which  controls  the  method  of  ray  energy 
summation.  A value  of  0 calls  for  incoherent  summation,  while  a value  of  1 
calls  for  semicoherent  summation.  The  program  also  permits  a value  of  2 to 
be  specified,  in  which  case  the  semicoherent  mode  of  summation  is  selected,  but 
without  the  attenuation  factor  which  cuts  down  the  amplitude  of  the  phase  in- 
terference oscillations  when  the  range  sampling  is  insufficient.  Use  of 
JCOH(IF)  = 2 is  not  recommended. 

Card  6 contains  the  source  depth  ZSO  in  the  first  field  and  the  NREC  re- 
ceiver depths  ZRO(IREC)  in  the  subsequent  fields.  These  depths  must  be  ex- 
pressed in  feet.  Neither  the  source  nor  the  receiver  may  be  located  within 
0.1  foot  of  the  surface  or  bottom. 

Card  7 contains  RMIN,  the  minimum  Tange,  DR,  the  range  increment,  and 
RMAX,  the  minimum  range,  all  in  yards.  The  number  of  range  points  is  not  per- 
mitted to  exceed  251.  If  the  specified  inputs  result  in  more  than  this  number, 
RMAX  is  automatically  truncated  to  the  maximum  allowable  value.  The  value  of 
the  parameter  KMAX,  which  represents  the  total  number  of  specified  range  points, 
is  evaluated  at  this  time. 

If  IBEAM  = 1,  the  beam  patterns  are  inserted  at  this  point  by  calling 
BEAM(O).  The  READ  statements  are  in  SUBROUTINE  BEAM(L).  The  patterns  aTe  read 
in  as  tables  of  pressure  ratio  (relative  to  the  beam  axis)  vs  ray  source  angle 
(degree),  beginning  at  -90  degrees  (downward)  and  extending  to  +90  degrees  (up- 
ward). NF  patterns  are  read  in,  one  for  each  frequency.  The  first  card,  8a, 
contains  the  integers  NBF(IF)  indicating  the  number  of  points  (up  to  a maximum 
of  100)  required  to  specify  each  of  the  patterns.  This  is  followed  by  card 
sets  8b  and  8c  for  the  first  frequency,  card  sets  8b  and  8c  for  the  second,  and 
so  on.  Card  set  8b  contains  the  grazing  angles  THETA(I)  and  card  set  8c  con- 
tains the  corresponding  beam  pressure  ratios  BM(I,IF).  The  beam  patterns  are 
applied  to  the  source  depth  ZSO  (referred  to  subsequently  as  the  true  source 
depth) . 

This  portion  of  PLRAY  concludes  with  a few  miscellaneous  calculations. 

The  volume  attenuation  coefficient  of  the  water  ALPHA(IF)  is  calculated  from 
Therp's  formula  (reference  (e))  for  each  of  the  frequencies.  SUBROUTINE  SURF 
is  called  to  compute  the  surface  reflection  loss  SLOSS(IF) . Except  for  name 
changes,  this  routine  is  identical  to  the  surface  loss  routine  of  reference 
(a).  Finally,  the  earth  curvature  correction  is  applied  to  the  source  depth, 
the  corrected  value  being  labeled  ZS. 

RECEIVER  DEPTH  LOOP 

The  remainder  of  the  program  is  enclosed  within  a receiver  depth  DO  loop. 
Since  only  one  receiver  depth  can  be  processed  at  a time,  it  is  necessary  to 
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carry  out  the  complete  set  of  calculations  for  one  depth  before  starting  on 
the  next.  At  the  beginning  of  the  loop  the  earth  curvature  correction  is  ap- 
plied to  the  current  receiver  depth,  the  corrected  value  being  labeled  ZR. 

Since  PLRAY  is  subject  to  large  errors  when  the  source  and  receiver  are  at  the 

same  depth,  ZS  and  ZR  are  compared  to  determine  whether  they  differ  by  at  least 

2 feet.  If  not,  ZS  is  automatically  shifted  upward  and  ZR  downward  by  a suf- 
ficient amount  to  assure  a 2-foot  separation.  The  next  step  is  the  setting  up 
of  a new  augmented  velocity  profile  by  inserting  the  source  and  receiver  depths 
into  the  old  one.  This  is  done  in  SUBROUTINE  INSERT. 

SUBROUTINE  INSERT 

The  depths  and  sound  speeds  of  the  original  profile  are  designated  ZB(IB) 
and  CB(IB),  the  index  being  IB.  The  principal  function  of  INSERT  is  to  set  up 

a new  table  of  ZA(IA)  and  CA(IA)  which  includes  the  source  and  receiver  depths. 

If  the  source  depth  happens  to  coincide  with  one  of  the  ZB(IB),  it  is  not  nec- 
essary to  add  a new  point.  Otherwise,  the  sound  speed  CS  is  calculated  at 
depth  ZS  from  the  appropriate  curvilinear  segment  parameters  and  this  point  is 
inserted  into  the  IA  table  at  the  proper  place.  The  indices  I A of  subsequent 
points  are  increased  by  1.  A similar  procedure  is  followed  for  the  receiver 
depth  ZR.  The  sound  speed  at  the  receiver  depth  is  designated  CR. 

Since  no  new  table  is  established  for  the  segment  parameters  CSP(IB) , 
GP(IB),  and  ZP(IB),  it  will  be  necessary  later  to  be  able  to  identify  the  val- 
ue of  IB  wh4  h corresponds  to  any  given  value  of  IA.  The  identification  is 
accomnl4  -etting  up  an  array  IDF(IA)  such  that 

i,  -vF(IA) 

INSERT  also  redefines  the  source  and  receiver  locations  for  use  in  the 
subsequent  ray  computations.  Let  ZAA  and  CAA  be  the  depth  and  sound  speed  of 
the  redefined  source  and  ZBB  and  CBB  the  depth  and  sound  speed  of  the  redefined 
receiver.  Then  CAA  is  defined  as  the  larger  of  CS  and  CR,  and  CBB  is  defined 
as  the  smaller.  Thus,  in  the  ray  computations  the  source  will  always  have  a 
larger  sound  speed  than  the  receiver.  To  avoid  confusion  in  terminology,  the 
original  source  and  receiver  will  henceforth  be  referred  to  as  the  true  source 
and  true  receiver.  Since  the  beam  pattern  is  always  applied  to  the  true  source, 
it  will  be  necessary  to  ascertain  whether  this  is  ZAA  or  ZBB.  The  identifica- 
tion parameter  is  SN.  It  is  assigned  a value  of  1 if  ZAA  is  the  true  source 
and  -1  if  ZBB  is  the  true  source. 

At  the  end  of  INSERT  is  a test  to  determine  whether  the  surface  duct  model 
is  to  be  implemented.  The  first  requirement  is  that  IDUCT  have  a value  other 
than  1,  that  is,  that  a duct  exist.  Secondly,  it  is  required  that  the  source 
be  located  within  the  duct  and  the  receiver  be  at  a depth  not  in  excess  of  1.8 
times  the  duct  depth,  or  vice  versa.  If  these  conditions  are  met,  the  param- 
eter IAD  is  set  to  the  value  of  the  index  IA  at  the  bottom  of  the  duct  (equiv- 
alent to  IDUCT  for  the  index  IB),  and  the  flag  ISD,  which  is  normally  zero,  is 
set  to  1. 

PROGRAM  PLRAY  (CONTINUED) 

Following  the  return  from  INSERT,  a value  of  1 or  2 is  assigned  to  the 
parameter  ICASE.  A value  of  1 applies  if  the  source  depth  ZAA  is  greater  than 
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the  receiver  depth  ZBB;  otherwise  ICASE  = 2.  This  parameter  is  required  be- 
cause of  the  difference  in  the  characteristics  of  the  multipath  ray  types,  de- 
pending upon  whether  the  source  is  above  or  below  the  receiver. 

SUBROUTINE  SECLIM 

At  this  point  SECLIM  is  called  for  the  purpose  of  setting  up  the  ray 
source  angle  sectors.  The  sectors  are  bounded  by  limiting  rays  which  separate 
rays  of  one  type  from  rays  of  another.  The  boundaries  of  the  sectors  are  de- 
fined in  terms  of  the  limiting  sound  speeds  CL(IL),  which  are  the  vertex  veloc- 
ities of  the  limiting  rays.  If  the  vertex  depth  of  the  limiting  ray  to  any 
sector  boundary  occurs  at  an  interior  maximum  point  of  the  velocity  profile, 
as  distinguished  from  the  top  or  bottom,  the  corresponding  value  of  CL(IL)  is 
flagged  with  a minus  sign. 

The  sound  speed  CL(1) , representing  the  inner  boundary  of  the  first  sector, 
is  the  vertex  velocity  of  the  shallowest  ray  which  reaches  the  receiver  depth. 
Since  the  sound  speed  at  the  source  is  by  definition  greater  than  (or  equal  to) 
the  sound  speed  at  the  receiver,  this  ray  is  normally  the  zero  degree  ray  and 
CL(IL)  is  equal  to  CAA.  However,  if  in  the  velocity  profile  there  exists  be- 
tween the  source  and  receiver  depths  a local  maximum  at  which  the  sound  speed 
exceeds  CAA,  CL(1)  is  set  to  the  negative  of  the  sound  speed  at  this  maximum. 

The  index  IA  of  the  depth  at  which  CL(1)  is  defined  is  identified  by  the  symbol 
ILT. 


Beginning  at  the  index  ILT,  a search  is  made  simultaneously  both  upward 
and  downward  along  the  velocity  profile  for  the  limiting  depths  which  define 
the  CL(IL)  at  the  sector  boundaries.  In  order  that  any  given  depth  ZA(IA)  may 
qualify  as  a suitable  limiting  depth,  it  is  necessary  first  of  all  that  the 
sound  speed  CA(IA)  exceed  the  absolute  value  of  the  previous  limiting  sound 
speed  CL(IL-l).  In  addition,  one  of  the  following  two  requirements  must  be 
met:  (1)  CA(IA)  must  occur  at  a maximum  of  the  profile,  that  is,  if  IA  = 1, 

CA(1)  must  exceed  CA(2) , if  IA  = NA,  CA(NA)  must  exceed  CA(NA-l) , and  otherwise 

CA(IA)  must  exceed  both  CA(IA-l)  and  CA(IA+1).  If  the  maximum  occurs  at  the 
surface  or  bottom,  CL(IL)  is  given  a positive  algebraic  sign;  otherwise  it  is 
given  a negative  sign.  (2)  If  the  above  condition  is  not  met,  tests  are  run 
on  the  sound  speed  gradient  G(IA),  which  is  calculated  at  the  beginning  of  the 
subroutine.  If  G(IA)  is  an  internal  maximum  or  minimum  value  (that  is,  not 
occurring  at  the  surface  or  bottom),  CL(IL)  is  set  equal  to  the  negative  of 
CA(IA) . 

Strictly  speaking,  the  ray  to  a segment  boundary  at  which  the  gradient  is 
an  extremum  (a  maximum  or  minimum)  is  not  a true  limiting  ray  and  does  not  form 
a boundary  between  families  of  rays  of  different  types.  However,  experience 

has  shown  that  the  presence  of  such  extrema  frequently  gives  rise  to  extrema  in 

the  ray  range,  i.e.,  to  caustics.  If  the  gradient  criterion  is  not  included, 
the  large  sector  formed  by  true  limiting  rays  sometimes  exhibits  multiple  caus- 
tics, which  cannot  be  handled  by  the  PLRAY  caustic  correction  procedure.  The 
gradient  criterion  has  been  inserted  to  divide  the  larger  sector  into  two 
smaller  ones,  thereby  alleviating  the  problem. 

After  all  the  values  of  CL(IL)  associated  with  the  velocity  profile  have 
been  determined,  a final  value  is  arbitrarily  assigned,  corresponding  to  a ray 
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which  leaves  the  source  at  30  degrees.  This  ray  forms  the  outer  boundary  of 
the  last  sector  in  which  rays  are  computed.  Beyond  30  degrees  the  ray  intensi- 
ties are  computed  from  an  approximate  formula  (appendix  F) . 

PROGRAM  PLRAY  (CONTINUED) 

Following  the  sector  determination,  PLRAY  initializes  several  variables 
prior  to  the  ray  computations.  The  H(K,IF)  array  and  the  parameters  ISURF  and 
IBOT  are  set  to  zero.  H(K,IF)  is  a large  array  in  which  the  various  ray,  caus- 
tic-correction and  surface  duct  intensities  are  stored  as  a function  of  range 
(K)  and  frequency  (IF). 

Before  entering  the  sector  loop,  PLRAY  tests  the  values  of  the  CL(IL)  to 
determine  which  sectors,  if  any,  contain  rays  which  propagate  within  a surface 
duct  and  assigns  to  the  variable  ILD  the  value  of  IL  at  the  outer  boundary  of 
the  highest  order  sector  of  this  type.  If  there  is  no  surface  duct,  the  value 
of  IAD  is  1. 

SECTOR  LOOP 

The  program  now  enters  a sector  DO  loop  (index  IL)  in  which  the  complete 
set  of  computations  are  carried  out  in  each  sector,  one  at  a time,  beginning 
with  the  innermost.  At  the  beginning  of  the  loop  the  angle  TH2  is  computed. 

TH1  and  TH2  normally  represent  the  angles  at  the  inner  and  outer  boundaries  of 
the  sector.  However,  on  the  first  transit  through  the  loop  (IL  = 1),  which 
serves  only  to  initialize  TH1,  the  value  assigned  to  TH2  is  the  angle  of  the 
inner  boundary  of  the  first  sector.  The  program  then  jumps  to  the  end  of  the 
loop  where  TH1  is  updated  by  setting  it  equal  to  TH2.  The  steady  state  is 
reached  on  the  second  transit  (IL  = 2),  where  TH1  and  TH2  have  been  assigned 
the  proper  values  corresponding  to  the  first  sector.  Thus,  each  sector  is 
identified  by  the  value  of  IL  at  its  outer  boundary. 

Immediately  following  the  computation  of  TH2,  the  index  IL  is  tested  to 
determine  whether  ray  computations  in  the  surface  duct  (if  any)  are  to  be  sup- 
pressed. If  IL  has  a value  less  than  or  equal  to  ILD  and  the  source  and  re- 
ceiver are  both  located  within  the  duct,  the  program  skips  to  the  end  of  the 
loop  and  updates  TH1  in  preparation  for  the  next  sector. 

Assuming  the  computations  in  the  sector  are  not  to  be  suppressed,  ISURF 
is  set  equal  to  1 if  the  rays  in  the  sector  strike  the  surface  and  IBOT  is  set 
to  1 if  the  rays  in  the  sector  strike  the  bottom. 

The  next  step  is  to  set  the  parameter  MJ.  MJ  is  a flag  used  to  identify 
a special  case  associated  with  the  first  sector  (IL  = 2).  Before  explaining 
the  function  of  MJ  it  will  be  instructive  to  discuss  the  way  in  which  the  var- 
ious multipath  ray  types  corresponding  to  any  given  ray  source  angle  are  nor- 
mally defined. 

MULTIPATH  RAY  TYPES 

The  bounding  angles  TH1  and  TH2  of  a sector  are  normally  defined  as  posi- 
tive numbers  and  the  individual  ray  source  angles  within  the  sector  are  normally 
defined  in  SUBROUTINE  DIVSEC  as  positive  numbers.  However,  the  multipath  ray 
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types  associated  with  a particular  source  angle  include  (where  applicable)  both 
up-going  and  down-going  rays.  The  ray  configurations  are  shown  schematically 
in  figures  la  and  lb  for  zero  order  arrivals  and  in  figures  2a  and  2b  for  first 
order  arrivals.  The  ray  configurations  for  the  second  and  higher  orders  are 
the  same  as  for  the  first,  the  only  difference  being  that  additional  ray  cycles 
are  added  to  the  four  existing  rays. 

In  the  subsequent  description  of  the  program,  the  different  multipath  ray 
types  will  be  identified  by  the  index  J.  As  may  be  seen  from  figures  la  and 
lb,  there  are  normally  two  multipath  ray  types  for  the  zero  order  arrivals. 

The  ray  SR1#  which  propagates  directly  from  the  source  to  the  receiver  with  no 
intervening  vertex,  is  identified  by  J = 1.  The  second  ray  SR?,  which  experi- 
ences one  upper  vertex,  is  identified  by  J = 2.  In  the  case  of  the  first  and 
subsequent  orders  there  are  normally  four  multipath  ray  types.  In  figures  2a 
and  2b  the  down-up  ray  SR^  is  identified  by  J = 1,  the  up-down-up  ray  SR2  by 
J = 2,  the  down-up-down  ray  SR3  by  J = 3,  and  the  up-down-up-down  ray  SR4  by 
J = 4. 

A more  realistic  portrayal  of  the  zero  order  rays,  exhibiting  the  refrac- 
tion resulting  from  the  profile  gradients,  is  shown  in  figures  3a,  3b,  and  3c. 
Figure  3a  is  a typical  example  of  ICASE  = 2.  It  will  be  observed  that  in  this 
case  the  rays  SR^  and  SR2  are  of  different  types;  considered  as  a function  of 
the  source  angle,  they  belong  to  different  families.  Figure  3b  is  an  example 
of  ICASE  = 1 in  which  there  is  a local  maximum  of  the  profile  between  the 
source  and  receiver  depths.  Here  again  the  rays  SR^  and  SR2  belong  to  separate 
families.  There  is  no  continuity  between  the  two  families  because  of  the  blank 
sector  formed  by  the  rays  which  are  trapped  in  the  surface  duct  and  never  reach 
the  receiver  depth. 

Consider  now  figure  3c.  In  this  case,  as  the  source  angle  is  reduced  to 
zero,  the  two  rays  SR}  and  SR2  become  coincident.  Thus,  the  up-going  rays  and 
down-going  rays  both  belong  to  one  family.  Therefore  in  this  case  the  bound- 
aries of  the  first  sector  are  redefined  so  as  to  include  both  negative  and  pos- 
itive source  angles.  The  lower  bounding  angle  TH1  is  set  equal  to  the  negative 
of  the  upper  bounding  angle  TH2,  and  there  is  now  only  one  multipath  ray  type 
(J  = 1)  in  the  zero  order  arrivals.  In  the  first  and  higher  order,  arrivals 
the  number  of  multipath  ray  types  is  reduced  from  4 to  2.  The  parameter  MJ  in- 
dicates the  number  of  multipath  ray  types  in  the  zero  order  arrivals  of  the 
first  sector,  two  normally,  and  one  in  the  above  special  case. 

SUBROUTINE  DIVSEC(IL) 

SUBROUTINE  DIVSEC  computes  the  source  angles  of  all  rays  to  be  computed  in 
the  sector.  To  avoid  computations  in  the  immediate  vicinity  of  limiting  rays, 
a small  dead  zone  of  width  DELO  = 0.00001  radian  is  inserted  at  each  of  the  two 
boundaries  of  the  sector,  so  that  the  source  angle  of  the  first  ray  computed  in 
the  sector  is  TH1  + DELO.  The  angular  spacing  DEL  between  successive  rays  is 
variable.  Beginning  at  a value  5.*DEL0,  it  is  successively  doubled  until  a 
limiting  value  DELM  is  reached.  The  value  of  DELM  is  0.01  radian  in  all  sec- 
tors but  the  last,  and  0.04  radian  in  the  last  sector.  The  reason  for  the 
larger  value  is  that  the  last  sector  consists  of  the  steeper,  less  profile-de- 
pendent bottom-bounce  rays,  for  which  the  range  vs  source  angle  curves  are 
well-behaved. 
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FIGURE  3a  - Multipath  Ray  Types,  Zero  Order  Arrivals,  ICASE  = 2 


FIGURE  3b  - Multipath  Ray  Types,  Zero  Order  Arrivals,  ICASE  = 1 


NADC-77296-30 


If  the  limiting  sound  speed  CL(IL)  at  the  outer  boundary  of  the  sector  is 
positive,  the  uniform  spacing  of  the  rays  is  continued  throughout  the  remainder 
of  the  sector.  However,  a negative  value  of  CL(IL)  signifies  that  the  outer 
boundary  is  formed  by  a limiting  ray  to  an  internal  maximum  of  the  profile. 

In  this  case,  as  the  outer  edge  of  the  sector  is  approached,  the  ray  range  ap- 
proaches infinity,  resulting  in  extremely  nonlinear  behavior.  Such  behavior 
requires  the  rays  to  be  more  closely  spaced.  In  the  subroutine  this  is  accom- 
plished by  applying  in  reverse  the  algorithm  used  in  the  vicinity  of  the  inner 
boundary. 

In  the  event  CL(IL)  is  negative,  the  ray  source  angles  in  the  sector  are 
defined  in  the  following  way.  Every  time  a source  angle  TH(I)  is  set  up,  a 
corresponding  temporary  angle  THT(I)  is  also  set  up  at  a corresponding  distance 
from  the  outer  boundary.  As  the  process  continues,  with  TH(I)  becoming  larger 
and  THT(I)  becoming  smaller,  eventually  the  two  angles  overlap.  As  soon  as  the 
overlap  occurs,  the  process  is  stopped  and  the  angles  THT  are  transferred  in  an 
appropriate  manner  to  the  TH  array.  The  total  number  of  rays  in  the  sector  is 
NTH. 


PROGRAM  PLRAY  (CONTINUED) 

The  basic  building  blocks  from  which  the  ranges  and  range  derivatives  are 
constructed  are  calculated  in  a sector  angle  (ITH)  DO  loop.  Before  calling 
SUBROUTINE  INCR  for  this  purpose,  PLRAY  computes  for  each  ray  in  the  sector 
the  vertex  velocity  CV  and  the  tangent  of  the  ray  source  angle  VA(ITH).  VA 
serves  as  the  independent  variable  for  the  range  derivatives,  in  lieu  of  the 
source  angle  itself. 

SUBROUTINE  INCR(ITH)  - COMPUTATION  OF  RANGE  AND  RANGE  DERIVATIVE  INCREMENTS 

In  the  first  portion  of  INCR  the  increments  of  the  range,  DX(IA),  and  the 
derivative  of  range  with  respect  to  VA,  DXP(IA),  are  computed  in  each  of  the 
applicable  layers  of  the  velocity  profile.  Then  these  individual  increments 
are  added  up  to  yield  the  three  basic  building  blocks,  namely,  (1)  the  range 
increment  RSR(ITH)  and  range  derivative  increment  RPSR(ITH)  corresponding  to 
the  portion  of  the  ray  path  between  the  source  and  receiver;  (2)  the  respective 
increments  RUP(ITH)  and  RPUP(ITH)  corresponding  to  the  portion  of  a ray  cycle 
which  lies  above  the  receiver  depth;  (3)  the  respective  increments  RCYC(ITH) 
and  RPCYC(ITH)  corresponding  to  a complete  ray  cycle.  The  formulas  from  which 
DX(IA)  and  DXP(IA)  are  computed  are  derived  in  appendix  G. 

Before  beginning  the  computations  it  is  necessary  to  determine  the  upper 
and  lower  vertex  depths.  To  determine  the  upper  vertex  depth,  a search  is  made 
through  all  the  layers  from  the  source  depth  to  the  surface  to  find  the  first 
layer  boundary  IA  at  which  the  sound  speed  CA(IA)  exceeds  the  vertex  velocity 
of  the  ray.  If  no  such  boundary  exists,  the  upper  vertex  occurs  at  the  sur- 
face. However,  if  the  boundary  is  found,  the  values  of  ZA(IA)  and  CA(IA)  are 
stored  temporarily  as  ZUT  and  CUT,  CA(IA)  is  replaced  with  the  vertex  velocity 
CV,  and  ZA(IA)  is  replaced  with  the  vertex  depth,  which  is  the  depth  on  the 
velocity  profile  corresponding  to  the  sound  speed  CV.  This  depth  is  computed 
by  solving  the  parabolic  formula  of  the  particular  profile  segment  in  which  CV 
lies. 
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A similar  procedure  is  followed  in  finding  the  lower  vertex  depth  by 
searching  downward  through  the  profile  layers  from  the  source  depth  to  the 
bottom.  If  a layer  boundary  IA  is  found  at  which  CA(IA)  exceeds  CV,  the  val- 
ues of  ZA(IA)  and  CA(IA)  are  stored  temporarily  as  ZDT  and  CDT  and  are  replaced 
with  the  lower  vertex  depth  and  the  vertex  velocity. 

The  increments  DX(IA)  and  DXP(IA)  are  computed  in  a DO  loop  in  which  the 
index  IA  runs  from  a minimum  value  IUP  determined  by  the  upper  vertex  depth  to 
a maximum  value  IDN  determined  by  the  lower  vertex  depth.  In  each  layer,  VI 
and  V2  represent  the  tangents  of  the  ray  angles  at  the  top  and  bottom  of  the 
layer  and  VV1  and  W2  are  their  respective  reciprocals.  If  the  ray  has  a re- 
fractive upper  vertex,  that  is,  if  the  ray  does  not  reach  the  surface,  VI  and 
W1  are  both  set  to  zero.  If  the  ray  has  a refractive  lower  vertex  (does  not 
reach  the  bottom),  V2  and  VV2  are  set  to  zero.  At  the  end  of  the  loop  the 
temporarily  stored  values  of  ZA  and  CA  are  replaced. 

The  increments  RSR(ITH),  RUP(ITH),  and  RCYC(ITH)  and  their  respective  de- 
rivatives are  then  computed  in  appropriate  DO  loops. 

Before  returning  to  PLRAY,  INCR  checks  IBOT  to  determine  whether  the  ray 
hits  the  bottom.  If  so,  the  ray  angle  at  the  bottom,  THBD  (degrees)  is  calcu- 
lated and  SUBROUTINE  BOTTOM  is  called  for  the  computation  of  the  bottom  loss 
BOTL(IF)  for  the  specified  frequencies.  These  values  are  then  stored  in  a 
two-dimensional  source  angle- frequency  array  BLOSS(ITH,IF) . 

SUBROUTINE  BOTTOM (THBD) 

This  subroutine  computes  the  bottom  loss  (db)  as  a function  of  grazing 
angle  THBD  and  frequency  F(IF).  The  bottom  loss  is  obtained  by  interpolating 
between  values  given  in  either  of  two  sets  of  tables,  either  the  internally 
stored  tables  based  on  the  FNWC  bottom  loss  curves,  or  the  external  tables  sup- 
plied by  the  user  in  the  input  data  deck.  If  the  bottom  parameter  IBP  is  as- 
signed a value  between  1 and  5,  the  FNWC  curves  are  used,  the  value  of  IBP  in- 
dicating which  of  the  five  FNWC  bottom  classes  is  to  be  selected.  If  IBP  = 6, 
the  bottom  loss  is  computed  from  the  externally  supplied  curves.  In  addition 
to  the  above,  two  other  options  are  available.  f IBP  = 0,  the  bottom  loss  is 
set  to  zero.  If  IBP  is  negative,  the  bottom  loss  is  arbitrarily  set  to  100  db. 
The  remainder  of  this  section  is  concerned  with  values  of  IBP  between  1 and  6. 

The  grazing  angles  and  corresponding  bottom  losses  of  the  tables  are 
DG(I,IBP,J)  and  DB(I,IBP,J),  where  the  index  I refers  to  the  point  number  along 
the  particular  bottom  loss  curve,  and  the  index  J identifies  the  frequency  to 
which  the  curve  applies.  The  actual  values  of  the  frequencies  of  the  bottom 
loss  curves  are  FR(J) . The  frequencies  at  which  the  propagation  loss  is  to  be 
computed  are  F(IF).  The  interpolations  are  carried  out  for  one  frequency  at  a 
time  (FF  = F(IF))  in  a frequency  DO  loop. 

In  the  general  case  a two-way  interpolation  is  required.  First,  the  two 
bottom  loss  curve  frequencies  FR(J)  and  FR(J+1)  are  selected  such  that  FF  lies 
between  them,  and  interpolations  with  respect  to  grazing  angle  are  carried  out 
to  obtain  the  two  bottom  loss  values  corresponding  to  THBD.  Then  an  interpola- 
tion with  respect  to  frequency  is  made  between  these  two  values  to  obtain  the 


I 

A 


28 


NADC- 77296- 30 


bottom  loss  corresponding  to  FF.  In  the  event  FF  coincides  with  one  of  the 
FR(J),  only  a single  interpolation  with  respect  to  grazing  angle  is  required. 

If  the  frequency  FF  is  less  than  FR(1) , it  is  assumed  that  the  bottom 
loss  curve  for  FR(1)  applies.  Similarly,  if  the  frequency  FF  is  larger  than 
FR(NFR) , it  is  assumed  that  the  bottom  loss  curve  for  FR(NFR)  applies.  As 
stated  previously,  the  frequencies  FR(J)  must  be  arranged  in  a monotonically 
increasing  order. 

PROGRAM  PLRAY  (CONTINUED) 

Having  calculated  the  basic  increments  of  range  and  range  derivative  for 
the  rays  of  the  sector,  the  program  is  now  ready  to  compute  the  ray  intensities 
and  their  multipath  sums.  The  computations  are  performed  first  for  the  zero 
order  arrivals,  following  which  a loop  is  entered  for  the  first  and  subsequent 
order  arrivals.  The  variable  indicating  the  order  of  the  arrivals  is  NCYC. 

The  number  of  multipath  ray  types  involved  is  NJ.  For  the  zero  order  arrivals 
NJ  = MJ;  for  the  higher  orders  NJ  = 2*MJ. 

ZERO  ORDER  ARRIVALS 

The  ranges  R(ITH,J)  and  range  derivatives  RP(ITH,J)  are  computed  from  the 
basic  increments  generated  by  INCR,  and  the  NJ  values  of  the  parameter  SIGN(J) 
are  determined.  SIGN(J)  is  the  algebraic  sign  of  the  angle  of  the  Jth  multi- 
path  ray  type  at  the  true  source.  When  a beam  pattern  is  specified,  these 
signs  are  needed  to  distinguish  between  up-going  and  down-going  rays  at  the 
true  source,  where  the  beam  pattern  is  applied.  The  remainder  of  the  computa- 
tions are  carried  out  in  RAYSUM  and  RCAUST.  Since  these  two  subroutines  are 
designed  to  process  all  arrivals,  regardless  of  the  order,  the  following  de- 
scription will  include  the  higher  orders  in  addition  to  the  zeroth. 

When  the  surface  duct  model  is  used  in  lieu  of  ray  theory  to  calculate  the 
propagation  within  a surface  duct,  it  has  been  found  necessary  to  suppress  not 
only  the  ducted  rays,  but  also  the  contributions  of  the  zero  order  arrivals  of 
all  rays.  Comparison  with  the  predictions  of  normal  mode  theory  shows  that  in- 
clusion of  the  contributions  of  the  steeper  rays  in  the  direct  propagation  zone 
leads  to  erroneously  high  intensities  at  short  ranges.  Although  originally 
overlooked  in  the  development  of  PLRAY,  this  result  is  reasonable,  since  the 
surface  duct  model  is  designed  to  cover  the  complete  propagation  in  the  duct, 
even  at  short  ranges.  The  computation  of  the  zero  order  arrivals  is  therefore 
bypassed  for  all  sectors. 

SUBROUTINE  RAYSUM(IL) 

Brief  Overview 

RAYSUM  processes  a single  arrival  order  of  the  rays  in  a single  sector. 
After  determining  the  minimum  and  maximum  ranges  covered  by  these  rays,  it 
searches  for  the  presence  of  caustics.  When  a caustic  is  found,  RAYSUM  calls 
RCAUST  for  the  computation  of  caustic-correction  intensities  over  an  interval 
of  ranges  on  either  side  of  the  caustic.  Except  for  a small  overlap  region 
where  both  caustic  and  ray  intensities  are  computed  and  compared  for  the  pur- 
pose of  providing  a smooth  transition  between  the  two  regions,  RAYSUM  bypasses 
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the  caustic-correction  region  and  computes  ray  intensities  over  the  remainder 
of  the  range  interval  covered  by  the  rays.  It  then  sums  the  intensities  of  the 
multipath  ray  types  at  each  range.  If  semicoherent  summation  is  requested, 
RAYSUM  computes  a coherence  factor  by  which  the  summed  intensity  at  each  range 
is  multiplied. 

Minimum  and  Maximum  Ranges 

RAYSUM  searches  through  the  ranges  of  all  the  rays  of  each  multipath  type 
in  the  sector  and  picks  out  the  minimum  and  maximum  range  RMNJ  and  RMXJ  of  each 
multipath  type  and  the  overall  minimum  and  maximum  range  RMN  and  RMX  among  all 
ray  types.  These  ranges  are  then  converted  to  values  of  the  range  index  K, 
namely,  KMNJ(J)  and  KMXJ(J)  for  each  individual  ray  type  and  KMN  and  KMX  for 
the  composite  of  all  ray  types.  Neither  KMX  nor  KMXJ(J)  may  exceed  the  limit- 
ing value  KMAX,  corresponding  to  the  maximum  specified  range. 

Search  for  Caustics 


The  search  is  carried  out  separately  for  each  multipath  ray  type  (each 
value  of  J) . A caustic  occurs  at  a point  where  the  range  vs  VA  (tangent  of  the 
source  angle)  curve  passes  through  a minimum  or  a maximum.  In  the  original 
version  of  the  program  the  procedure  consisted  simply  of  a search  for  one  or 
the  other.  Experience  has  shown,  however,  that  sometimes  these  curves  are  not 
simply  parabola-like,  but  may  exhibit  more  complicated  behavior.  Although  it 
has  proved  impossible  to  take  care  of  all  eventualities,  an  algorithm  has  been 
inserted  to  deal  with  one  of  the  more  common  complications. 

Figure  4 shows  a range  curve  exhibiting  both  a maximum  and  a minimum. 
Although  it  may  be  possible  to  include  two  caustic  corrections  in  one  interval, 
the  complications  in  the  program  become  severe.  Fortunately,  in  most  cases 
examined,  one  of  the  two  extrema  has  been  found  to  be  a relatively  small  dip 
or  rise  at  one  end  of  the  curve.  As  currently  written,  the  program  searches 
for  the  more  important  of  the  two  extrema  and  ignores  the  other.  The  price 
paid  for  neglecting  the  other  extremum  is  an  occasional  caustic  spike  in  the 
resultant  propagation  loss  curve. 

Referring  to  figure  4,  RMXJ  and  RMNJ  are  the  largest  and  smallest  ranges 
in  the  sector.  RU  and  RD  are  the  larger  and  smaller  of  the  ranges  of  the  first 
and  last  rays  computed  in  the  sector.  The  significance  of  the  range  differences 
DRUP  and  DRDN  is  apparent  in  the  figure.  If  there  is  no  extremum  in  the  curve, 
DRUP  and  DRDN  are  both  zero.  If  there  is  only  a maximum,  DRDN  is  zero.  If 
there  is  only  a minimum,  DRUP  is  zero.  When  both  extrema  occur,  DRUP  and  DRDN 
are  compared  and  the  extremum  associated  with  the  smaller  of  the  two  is  ignored. 
In  the  example  shown  in  figure  4 the  minimum  is  ignored  and  the  correction  is 
applied  only  to  the  maximum. 

The  parameter  ICAUST(J)  is  a flag  indicating  whether  the  Jth  multipath  ray 
type  in  the  sector  experiences  a caustic  and,  if  so,  whether  it  is  a minimum- 
or  maximum- range  caustic.  A value  of  0 signifies  no  caustic,  1 a minimum-range 
caustic,  and  -1  a maximum-range  caustic.  If  a caustic  is  found,  the  value  of 
the  index  ITH,  such  that  the  ray  to  the  caustic  lies  between  ITH-1  and  ITH,  is 
stored  as  ICTH(J).  When  a caustic  is  present,  SUBROUTINE  RCAUST  is  called. 

The  parameter  KT1I  in  the  call  is  equivalent  to  ICTH(J). 
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SUBROUTINE  RCAUST ( KTH , J ) 

A derivation  of  the  caustic  correction  formulas  will  be  found  in  appendix 
H.  The  key  parameter  determining  the  intensity  at  a caustic  is  the  second  de- 
rivative of  the  range,  RPP.  In  this  program  the  derivative  is  taken  with  re- 
spect to  V A.  A tentative  value  of  RPP  is  computed  from  the  average  slope  of 
the  first  derivative  between  the  points  KTH-1  and  KTH  on  either  side  of  the 
caustic,  the  assumption  being  that  in  the  vicinity  of  the  caustic  the  range 
curve  is  parabolic  and  the  first  derivative  linear.  The  value  VAC  of  the  tan- 
gent function  VA  at  the  caustic  is  computed  from  the  first  and  second  range 
derivatives  on  the  same  assumption.  Then  two  values  of  the  range  to  the  caus- 
tic RC  are  computed  on  the  assumption  of  parabolic  behavior,  one  based  on  the 
point  at  KTH-1  and  the  other  at  KTH;  the  average  of  the  two  is  the  range  RC. 

It  has  been  observed  in  a number  of  instances  that  the  actual  behavior  of 

the  range  derivative  curves  deviates  rather  badly  from  the  straight-line  as- 
sumption upon  which  the  caustic  correction  theory  is  based.  The  difficulty 
tends  to  show  up  chiefly  in  an  abnormally  high  value  of  the  second  derivative 

RPP.  Experience  has  shown  that  more  reasonable  results  can  usually  be  ob- 

tained by  making  alternate  calculations  of  RPP  from  the  range  derivatives  cor- 
responding to  adjacent  pairs  of  rays  on  either  side  of  the  caustic,  namely, 
KTH-2,  KTH-1, and  KTH,  KTH+1,  and  selecting  the  numerically  smallest  of  the 
three  values. 

The  first  portion  of  RCAUST  continues  with  the  computation  of  numerous 
frequency- independent  quantities  required  for  the  caustic  correction.  If  a 
beam  pattern  has  been  specified,  it  is  necessary  to  compute  VBM,  the  tangent 
of  the  angle  of  the  ray  to  the  caustic,  evaluated  at  the  depth  of  the  true 
source.  The  flag  SN  determines  whether  the  true  source  is  the  ray  source  or 
the  ray  receiver  and  SIGN(J)  determines  whether  this  ray  leaves  the  true  source 
in  an  upward  or  a downward  direction.  The  beam  pattern  factors  (pressure  ra- 
tios) BMFAC(IF)  are  computed  in  SUBROUTINE  BEAM,  which  will  be  described  later. 
The  intensities  throughout  the  caustic-correction  range  interval  are  multiplied 
by  the  square  of  BMFAC(IF),  computed  from  the  value  of  VBM  at  the  caustic. 

Frequency  Loop 

The  remainder  of  RCAUST  consists  of  a frequency  DO  loop  in  which  the  in- 
tensities on  either  side  of  the  caustic  are  computed,  one  frequency  at  a time. 

Range  Limits  of  the  Caustic  Correction 

The  caustic  correction  intensity  is  proportional  to  the  square  of  the  Airy 
function  Ai(x),  the  constant  of  proportionality  being  COEF2,  a coefficient 
which  varies  inversely  with  the  2/3  power  of  RPP.  The  argument  x of  the  Airy 
function  is  proportional  to  the  horizontal  distance  from  the  caustic  at  the  re- 
ceiver depth.  The  constant  of  proportionality,  C0EF1,  varies  inversely  with 
the  1/3  power  of  RPP.  On  the  shadow  side  of  the  caustic,  where  x is  positive, 
Ai(x)  decays  exponentially.  An  arbitrary  limit  of  3.5  is  imposed  on  x as  the 
point  where  the  intensity  has  decayed  to  a negligible  value.  The  range  cor- 
responding to  this  limiting  value  of  x is  RC1.  For  a minimum  range  caustic 
RC1  is  less  than  RC.  For  a maximum  range  caustic  it  is  greater. 
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On  the  illuminated  side  of  the  caustic,  where  x is  negative,  the  Airy 
function  is  oscillatory.  On  this  side,  two  limiting  values  are  assigned  to  x. 
The  first  limit,  x = -1.77,  corresponds  (in  an  ideal  case  of  a parabolic  range 
curve)  to  a point  beyond  the  first  maximum  where  the  intensity  is  equal  to  the 
incoherent  sum  of  the  two  rays  which  form  the  caustic.  Ideally,  termination 
of  the  caustic  corrections  at  this  limit  assures  continuity  with  the  adjacent 
ray  intensities.  In  working  with  actual  velocity  profiles,  however,  conditions 
are  never  ideal,  and  a discontinuity  is  usually  found,  with  the  caustic-correc- 
tion intensity  too  high  relative  to  the  ray  intensity.  To  smooth  the  discon- 
tinuity, the  caustic-correction  computations  are  extended  out  to  a second  lim- 
it, x = -2.2,  which  is  close  to  the  first  zero  of  the  Airy  function.  The  sum 
of  the  two  ray  intensities  is  also  computed  in  the  overlap  interval  from  -1.77 
to  *2.2  and  continuity  is  assured  by  selecting  the  larger  of  the  two  intensi- 
ties at  each  range  point  in  the  interval.  The  range  corresponding  to  the  first 
limit  -1.77  is  RC2  and  the  range  corresponding  to  the  second  limit  -2.2  is  RC3. 
For  a minimum  range  caustic  RC2  and  RC3  are  larger  than  RC,  and  for  a maximum 
range  caustic  they  are  smaller. 

Range  indices  Kl,  K2,  and  K3  corresponding  to  RC1,  RC2 , and  RC3  respec- 
tively are  established  such  that  the  set  of  values  of  K from  Kl  to  K2  inclusive 
contains  all  the  range  points  between  RC1  and  RC2,  and  the  set  of  values  from 
(but  not  including)  K2  to  and  including  K3  contains  all  the  range  points  be- 
tween RC2  and  RC3.  As  a minor  refinement  to  the  above,  the  value  of  K3  is  al- 
tered, if  necessary,  so  that  the  difference  between  K2  and  K3  does  not  exceed 
15. 

Intensities 


The  caustic-correction  intensity,  identified  by  the  symbol  TRM,  consists 
of  two  factors,  HK  and  FAC.  HK  is  associated  with  the  caustic  itself  and  cor- 
responds to  the  spreading  loss  factor  in  the  ray  intensity.  FAC  is  an  atten- 
uation factor  which  includes  the  volume  attenuation  of  the  water,  the  bottom 
and  surface  losses,  and  the  beam  deviation  loss.  The  beam  loss  is  expressed 
in  the  form  of  an  intensity  ratio,  which  is  the  square  of  the  pressure  ratio 
BMFAC(IF) . The  intensities  computed  in  the  interval  between  Kl  and  K2  are 
added  directly  to  the  H(K,IF)  array.  The  intensities  computed  in  the  overlap 
interval  from  K2  to  K3  are  stored  temporarily  in  an  array  HC(KC,J,IF)  for  later 
comparison  with  ray  intensities  in  the  same  interval.  KC  is  a relative  range 
index  which  begins  at  K2+1  for  a minimum  range  caustic  and  at  K3  for  a maximum 
range  caustic. 

Since  the  limits  K2  and  K3  are  needed  in  RAYSUM  to  define  the  interval 
over  which  the  intensity  comparisons  are  to  be  made,  they  are  transmitted  in 
COMMON  through  the  variables  KA(J,IF)  respectively. 

SUBROUTINE  RAYSUM ( I L)  (CONTINUED) 

Parameters  for  Simplified  Computations  at  Steep  Angles 

Upon  the  return  from  RCAUST  a check  is  made  to  determine  whether  the  cur- 
rent sector  is  the  final  one  (IL  = NL) . If  so,  the  parameters  A0(J)  and  A1 ( J) 
are  evaluated  for  each  value  of  J.  These  parameters  are  required  for  the  steep 
angle  computation «. 
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Range  Loop 

The  remainder  of  RAYSUM  is  contained  within  a range  loop  in  which  the  in- 
dex K runs  from  a lower  limit  K1  to  the  already  determined  upper  limit  KMX. 

In  all  sectors  but  the  last,  K1  has  the  normal  lue  KMN.  The  last  sector  con- 
tains the  steepest  rays,  extending  out  to  30  J<  rees.  The  minimum  value  of  K 
for  these  rays  is  of  course  KMN.  However,  pro  j.sion  is  made  in  the  final  sec- 
tor for  the  computation  by  a simplified  approximate  method  (appendix  F)  of  the 
intensities  at  shorter  ranges,  where  K runs  from  1 to  KMN-1.  For  this  reason 
K1  is  assigned  a value  of  1. 

The  Index  JJ 


The  presence  of  a minimum  or  maximum  in  the  range  curve  leads  to  an  addi- 
tional complication  besides  the  problem  of  caustics.  The  complication  arises 
from  the  existence  of  two  branches  on  the  range  curve,  so  that  there  are  now 
two  rays  with  different  source  angles  propagating  to  the  same  range  point. 

This  situation  is  illustrated  in  figure  5,  which  shows  a range  curve  containing 
a minimum.  The  shaded  area  surrounding  the  minimum  point  is  the  region  uncon- 
ditionally covered  by  the  caustic  correction.  Above  this  is  a second  shaded 
area  in  which  both  caustic  and  ray  intensities  are  computed.  The  unshaded  area 
at  the  top  is  the  region  covered  by  ray  computations  only.  For  the  purpose  of 
the  ray  intensity  computations  the  two  branches  of  the  curve  are  identified  by 
different  values  of  the  index  JJ.  On  the  left  branch,  which  has  a negative 
slope,  JJ  is  assigned  a value  1,  while  on  the  right  branch,  which  has  a posi- 
tive slope,  the  value  is  2. 

The  presence  of  the  two  branches  requires  the  setting  up  of  three  new  ar- 
rays, VAKK(J.JJ)  for  the  tangent  function,  RPKK(J,JJ)  for  the  range  derivative, 
and  HTRM(J,JJ,IF)  for  the  ray  intensity.  The  VAKK  and  HTRM  arrays  are  set  to 
zero  for  each  value  of  K at  the  beginning  of  the  range  loop. 

Ray  Intensities 

The  ray  intensity  computations  are  performed  within  a DO  loop  over  the  ray 
type  index  J.  At  each  specified  range  point  the  intensities  are  computed  se- 
quentially for  each  multipath  ray  type.  Before  continuing  with  a description 
of  the  program  details,  it  will  be  useful  to  present  a brief  overview  of  the 
manner  in  which  the  ray  intensities  are  computed. 

The  formula  for  the  ray  intensity  consists  of  the  product  of  a spreading 
loss  factor  and  an  attenuation  factor  which  includes  the  losses  due  to  surface 
and  bottom  reflections,  volume  attenuation  of  the  water,  and  (where  applicable) 
the  beam  pattern.  The  formula  for  the  spreading  loss  factor  is  discussed  in 
appendix  I. 

The  range-dependent  variables  which  enter  into  the  formula  for  the  inten- 
sity include  VAK,  the  tangent  of  the  source  angle  of  the  ray  which  propagates 
to  the  range  K,  VBK,  the  tangent  of  the  angle  at  the  receiver,  RPK,  the  range 
derivative  (actually  RPKA,  its  magnitude),  CVK,  the  ray  vertex  velocity  (actu- 
ally CVK2,  its  square),  and  BLK,  the  bottom  loss.  The  primary  variables  are 
VAK,  RPK,  and  BLK;  CVK  and  VBK  are  derived  from  VAK.  The  letter  K appearing 
in  these  variable  names  signifies  that  they  pertain  to  the  ray  which  propagates 
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to  the  range  point  with  index  K,  identified  as  RK.  Since  the  ray  has  not  ac- 
tually been  computed,  it  is  necessary  to  evaluate  the  primary  variables  from 
the  information  at  hand. 

In  the  normal  method  of  computing  intensities  these  variables  are  evalu- 
ated by  interpolating  between  adjacent  rays  of  the  family  computed  in  the  cur- 
rent sector.  These  rays  are  identified  by  the  index  ITH.  The  procedure  is  to 
search  through  the  family  until  a value  of  ITH  is  found  such  that  the  range 
point  RK  lies  between  the  adjacent  ranges  R(ITH-1,J)  and  R(ITH,J).  VAK,  RPK, 
and  BLK  are  then  obtained  by  linear  interpolation  in  the  VA,  RP,  and  BLOSS  ar- 
rays. In  those  cases  where  no  caustic  is  present  and  the  R vs  VA  curve  has  no 
extremum,  the  interpolation  is  accomplished  by  a single  search  through  all  the 
rays  of  the  sector.  When  a caustic  is  present  the  curve  has  two  branches,  and 
it  is  necessary  to  investigate  each  branch  separately.  To  cover  all  cases  in 
a single  loop,  the  minimum  and  maximum  values  of  ITH  between  which  the  search 
is  conducted  are  denoted  by  ITHMN  and  ITHMX.  The  number  of  interpolations  is 
controlled  by  the  variable  NLOOP. 

The  normal  method  of  computing  intensities  applies  to  all  range  points 
covered  by  the  family  of  rays  in  each  sector,  the  range  index  K running  from 
KMNJ(J)  to  KMXJ(J).  In  the  final  sector,  the  simplified  approximation  is  em- 
ployed to  compute  VAK  and  RPK  for  propagation  at  angles  steeper  than  30  de- 
grees, where  no  rays  are  actually  computed.  The  remaining  variables  required 
for  the  intensity  are  then  derived  from  VAK  and  RPK. 

Let  us  now  return  to  the  program.  At  the  beginning  of  the  J loop  the 
variable  KTYP  is  set  to  zero.  KTYP  is  a flag  which  indicates  whether  the  in- 
tensity is  to  be  computed  in  the  normal  manner  (KTYP  = 0)  or  from  the  simpli- 
fied formula  for  steep  angles  (KTYP  = 1) . ILOOP  and  NLOOP  are  set  to  1 and 
ITHMN  and  ITHMX  are  set  to  1 and  NTH,  respectively.  A test  is  then  made  to 
determine  whether  the  special  steep-angle  method  is  to  be  employed.  If  so, 

KTYP  is  set  to  1 and  VA  and  RPK  are  computed  from  the  simplified  formulas 
(appendix  F) . Since  the  bottom  loss  has  not  been  precomputed  in  this  region, 
it  is  necessary  to  compute  the  bottom  grazing  angle  and  call  the  bottom  sub- 
routine . 

If  the  intensity  is  to  be  computed  in  the  normal  manner,  ICAUST(J)  is 
tested  to  determine  whether  a caustic  is  present.  If  not,  the  program  jumps 
to  the  JTH  DO  loop  in  which  the  search  procedure  is  carried  out. 

If  a maximum  range  caustic  is  present  (ICAUST(J)  = -1),  a test  is  made  to 
determine  whether  the  range  RK  lies  within  the  caustic  correction  region  for 
the  highest  frequency  being  processed.  If  so,  no  ray  intensity  computation  is 
to  be  computed  for  any  of  the  frequencies.  If  not,  tests  are  made  to  determine 
whether  the  range  RK  lies  within  either  or  both  of  the  two  branches  of  the 
range-VA  curve.  The  branches  will  be  labeled  left  and  right  on  the  assumption 
that  the  index  ITH  increases  from  1 at  the  extreme  left  to  NTH  at  the  extreme 
right. 

If  RK  lies  only  on  the  left  branch,  ITHMX  is  set  to  ICTH(J) , the  index  of 
the  ray  adjacent  to  the  caustic.  ILOOP  and  NLOOP  retain  their  original  values 
of  1.  If  RK  lies  only  on  the  right  branch,  ITHMN  is  set  to  ICTH(J)  and  ILOOP 
and  NLOOP  are  both  changed  to  2 (indicating  the  second  branch).  If  RK  lies  on 
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both  branches,  NLOOP  is  changed  to  2 while  I LOOP  remains  at  1,  thereby  signify- 
ing that  there  will  be  two  searches  and  interpolations.  ITHMX  is  set  to  ICTH(J) 
in  preparation  for  the  first  search,  which  occurs  on  the  left  branch. 

When  a minimum  range  caustic  is  present,  ICAUST(J)  being  positive,  a sim- 
ilar procedure  is  followed,  with  appropriate  interchanging  of  the  "greater 
than"  and  "less  than"  symbols. 

The  search  procedure  is  carried  out  in  the  JTH  DO  loop.  When  no  caustic 
is  present  the  rays  are  examined  in  increasing  order  from  1 to  NTH.  When  a 
caustic  is  present,  the  search  begins  at  ICTH(J)  and  proceeds  outward  toward 
the  edges.  After  the  proper  value  of  ITH  has  been  determined,  the  interpola- 
tions for  VAK  and  RPK  are  performed. 

The  procedure  from  this  point  on  applies  to  both  the  normal  case  and  the 
special  steep-angle  case.  The  algebraic  sign  of  the  range  derivative  RPK  is 
tested  and  the  appropriate  value  of  the  parameter  JJ  is  assigned,  1 if  RPK  is 
negative  and  2 if  positive.  The  value  of  VAK  is  stored  in  the  appropriate  bin 
of  the  VAKK(J,JJ)  array  and  the  absolute  value  RPKA  is  computed  and  stored  in 
the  RPKK(J,JJ)  array.  Finally,  the  spreading  loss  factor  HK  of  the  ray  inten- 
sity is  computed. 

It  now  remains  to  compute  the  frequency-dependent  factors.  First,  if  a 
beam  has  been  specified,  the  beam  pressure  ratios  BMFAC(IF)  are  obtained  by 
calling  BEAM(l).  The  angular  input  from  which  these  functions  are  computed  is 
VBM,  the  tangent  of  the  ray  angle  at  the  true  source.  VBM  is  obtained  from 
VAK  if  the  true  source  is  located  at  the  ray  source  (SN  = +1.)  and  from  VBK  if 
the  true  source  is  located  at  the  ray  receiver  (SN  = -1.).  SIGN(J)  indicates 
whether  the  ray  is  up-going  (+)  or  down-going  (-).  The  values  of  BMFAC(IF) 
are  stored  in  the  BMF(J,JJ,IF)  array  for  future  use  in  the  coherence  computa- 
tions. If  no  beam  has  been  specified,  these  steps  are  bypassed. 

The  remainder  of  the  intensity  calculations  are  performed  in  a frequency 
DO  loop.  Before  proceeding  further  it  is  necessary  to  ascertain  whether  the 
intensity  corresponding  to  this  particular  range  and  frequency  has  already  been 
calculated  in  SUBROUTINE  RCAUST.  A partial  exclusion  has  already  been  made  by 
comparing  the  range  index  K with  the  caustic  correction  limit  KA(J.IFMAX)  for 
the  highest  frequency.  At  lower  frequencies  the  extent  in  range  of  the  caus- 
tic-correction region  increases.  K is  now  compared  with  KA(J,IF)  for  the  cur- 
rent frequency.  If  the  test  shows  that  the  present  range- frequency  point  has 
not  been  covered  in  RCAUST,  the  attenuation  factor  FAC  is  calculated.  FAC  is 
composed  of  the  beam  intensity  ratio,  which  is  the  square  of  BMFAC(IF) , and  the 
antilogs  of  the  bottom  loss  BLK,  surface  loss  SLK,  and  the  volume  attenuation 
loss.  If  the  special  steep-angle  method  has  been  used  (KTYP  = 1),  the  bottom 
loss  has  already  been  computed.  Otherwise  (KTYP  = 0)  it  is  necessary  to  in- 
terpolate between  the  appropriate  pair  of  rays  in  the  BLOSS  array.  The  result- 
ing intensity  is  stored  in  the  HTRM(J,JJ, IF)  array. 

At  the  conclusion  of  the  intensity  calculations  a test  of  I LOOP  is  made 
to  ascertain  whether  the  procedure  must  be  repeated  for  the  other  branch  (if 
any)  of  the  range-VA  curve. 


37 


NADC-77296-30 


« 


There  remains  yet  one  more  item  of  business  before  completion  of  this 
section  of  P.AYSUM  - comparison  of  the  ray  intensities  with  the  caustic-corTec- 
tion  intensities  in  the  overlap  region.  The  comparison  is  made  separately  for 
each  frequency.  At  the  outset  a test  is  made  to  ascertain  whether  the  current 
range- frequency  point  lies  within  the  transition  region.  If  so,  the  intensi- 
ties of  the  rays  on  the  two  branches,  HTRM(J,1,IF)  and  HTRM(J,2,1F)  are  added 
to  yield  TRM  and  TRM  is  compared  with  the  caustic  intensity  HC(KC,J,IF).  If 
HC(KC,J,IF)  is  the  larger,  it  is  added  to  the  H(K,IF)  array  and  both  HTRM(J,1,IF) 
and  HTRM(J,2,IF)  are  set  to  zero.  If  TRM  is  the  larger,  nothing  is  added  to 
the  H(K,IF)  array  at  this  point.  Instead,  the  two  HTRM  values  are  retained  for 
use  in  the  final  summation. 

Ray  Summation 

The  remainder  of  SUBROUTINE  RAYSUM  is  devoted  to  the  summation  of  the  mul- 
tipath ray  intensities  and  computation  of  coherence  factors.  Since  it  is  as- 
sumed that  there  is  no  coherence  between  rays  belonging  to  different  branches 
of  the  range- VA  curve,  rays  corresponding  to  the  two  different  values  of  JJ  may 
be  processed  separately.  For  this  reason  the  ray  summation  and  coherence  com- 
putations are  made  in  a DO  loop  over  JJ  which  extends  to  the  end  of  the  sub- 
routine. Nested  inside  the  JJ  loop  is  a frequency  loop  which  also  extends  to 
the  end. 

At  the  outset  the  coherence  parameter  ICOH  and  the  variable  TRM  are  ini- 
tialized to  zero.  A zero  value  of  ICOH  signifies  that  no  coherence  exists  at 
either  end  of  the  rays.  It  will  later  be  set  to  1 if  coherence  occurs  only  at 
the  source  end,  to  2 if  only  at  the  receiver  end,  and  to  3 if  at  both  ends. 

The  variable  TRM  will  denote  the  resultant  intensity  of  the  current  arrival 
order,  ultimately  to  be  added  to  the  appropriate  bin  of  the  H(K,IF)  array. 

The  intensities  of  the  individual  rays  are  then  summed  in  a loop  over  the 
multipath  index  J to  yield  the  resultant  incoherent  intensity  HK  of  the  arrival 
order.  In  the  same  loop  are  computed  the  number  of  contributing  rays,  ISUM, 
and  the  code  parameter  IND,  which  indicates  exactly  which  rays  are  contribut- 
ing. The  J values  of  the  contributing  rays  are  arranged  in  IND  in  increasing 
order  from  left  to  right.  It  should  be  noted  that  it  is  not  necessarily  true 
that  all  of  the  NJ  rays  contribute  to  the  multipath  sum.  It  is  possible  that 
the  receiver  location  may  be  reached  by  some  of  the  rays  but  lie  in  a shadow 
zone  of  others.  It  is  also  possible  that  one  or  more  of  the  rays  may  lie  in  a 
range  interval  covered  by  a caustic  correction  and  therefore  be  excluded  from 
the  ray  intensity  computation. 

At  this  point  if  ISUM  is  zero,  the  program  jumps  to  the  end  of  the  fre- 
quency loop,  since  there  are  no  contributions  to  the  intensity.  If  incoherent 
summation  has  been  requested  (JCOH(IF)  = 1)  or  if  the  rays  do  not  reach  the 
surface  (ISURF  = 0)  or  if  there  is  only  one  contributing  ray  (ISUM  = 1) , no 
further  computations  are  necessary  and  the  value  of  TRM  (zero  in  this  case)  is 
augmented  by  HK  in  preparation  for  insertion  into  the  final  H(K,KF)  array.  If 
ISUM  has  a value  of  2,  3,  or  4,  the  program  proceeds  to  ascertain  whether  co- 
herence occurs  and,  if  so,  to  compute  the  coherence  factors. 
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Coherence  Computations 

A discussion  of  the  background  of  the  coherence  feature  of  PLRAY  and  a 
derivation  of  the  formulas  for  the  coherence  factors  is  given  in  appendix  J. 
Coherence  factors  are  computed  either  for  pairs  or  rays  or  for  all  four  rays. 

In  the  case  of  two  rays  (ISUM  = 2) , a test  is  made  to  determine  whether  ICOH 
is  0,  1,  or  2.  If  it  is  zero  there  is  no  coherence  and  the  resultant  inten- 
sity is  simply  the  incoherent  sum.  If  ICOH  = 1,  the  coherence  factor  is  com- 
puted at  the  ray  source  end.  If  it  is  2,  the  coherence  factor  is  computed  at 
the  ray  receiver  end.  If  ISUM  = 3,  one  pair  of  rays  is  selected  and  processed 
as  in  the  two-ray  case,  and  the  third  is  added  incoherently.  If  ISUM  = 4,  ICOH 
may  have  a value  of  0,  1,  2,  or  3.  If  ICOH  = 0,  the  resultant  intensity  is 
just  the  incoherent  sum.  If  1 or  2,  the  four  rays  are  split  into  two  pairs 
and  each  pair  is  processed  as  in  the  two-ray  case.  If  ICOH  = 3,  the  coherence 
factors  are  computed  for  the  combined  interference  of  all  four  rays.  When  two 
pairs  of  rays  are  to  be  processed,  the  loop  in  which  the  processing  is  done 
must  be  traversed  twice.  The  parameter  which  controls  the  number  of  transits 
through  the  loop  is  NLOOP. 

Consider  now  the  case  of  two  rays.  If  ISUM  = 2,  a check  is  made  to  see 
if  the  contributing  rays  are  1 and  4 or  2 and  3,  in  which  case  coherence  is 
assumed  not  to  occur.  Otherwise,  the  average  value  of  VAK  for  the  two  rays  is 
computed,  and  from  this  the  average  vertex  velocity  CVK  and  the  average  tangent 
of  the  ray  angle  at  the  surface,  VO.  It  is  necessary  at  this  point  to  deter- 
mine whether  the  coherence  is  to  be  considered  at  the  source  or  receiver  end. 
The  source  end  is  selected  for  zero  order  arrivals  if  ZAA  <ZBB  and  for  first 

and  higher  order  arrivals  if  the  contributing  rays  are  1 and  2 or  3 and  4. 

Otherwise  the  receiver  end  is  selected.  The  horizontal  separation  of  the  rays 
ZRCOH  is  computed  from  equation  (J-l)  of  appendix  J and  computed  with  the 
threshold  value  RCOH.  If  the  threshold  is  exceeded,  ICOH  remains  zero  and  the 
incoherent  sum  applies.  Otherwise  ICOH  is  set  to  1 if  the  coherence  is  at  the 
source  end  and  to  ICOH  + 2 if  the  coherence  is  at  the  receiver  end.  NLOOP  is 
set  to  1. 

The  phase  angle  is  now  computed  from  equation  (J-3) , using  ZAA  or  ZBB  for 

the  depth,  according  to  the  value  of  ICOH,  and  the  coefficient 

FRACTS  - 2gagb/(8a2.gb2) 

is  also  computed,  where  ga  and  g^  are  the  beam  pattern  functions  (pressure 
ratios)  of  the  two  rays,  their  FORTRAN  equivalents  being  BMF1  and  BMF2.  If  no 
beam  has  been  specified,  FRACTS  is  given  the  value  of  1.  If  "full"  coherence 
(JCOH  = 2)  has  been  requested,  the  coherence  factor  COFACS  (second  factor  in 
equation  (J-5))  is  computed  and  multiplied  by  the  incoherent  sum  HK  to  yield 
the  resultant  coherent  intensity.  However,  the  "full"  coherence  option  is  not 
recommended  for  normal  use.  If  the  semicoherent  option  (JCOH  = 0)  has  been 
specified,  it  is  necessary  to  consider  the  sampling  problem.  The  attenuation 
factor  FRACTB  (F  in  appendix  J)  is  computed  from  equation  (J-15)  and  compared 
with  FRACTA,  the  absolute  value  of  FRACTS.  If  FRACTB  is  smaller,  FRACTS  is 
replaced  by  FRACTB  with  the  algebraic  sign  of  FRACTS  appended. 


NADC-77296-30 


Consider  next  the  case  of  three  rays.  If  ISUM  = 3,  one  of  the  three  rays 
is  treated  incoherently  and  the  other  two  rays  are  handled  as  in  the  two-ray 
case.  The  rules  for  selecting  the  rays  are  given  in  table  J-I  of  appendix  J. 
Prior  to  processing  the  pair  of  rays  the  intensity  of  the  single  ray  must  be 
subtracted  from  the  three-ray  sum  and  included  instead  in  the  variable  TRM. 

Consider  finally  the  case  of  four  rays.  To  determine  the  value  of  ICOH, 
the  average  of  all  four  values  of  VAK  is  computed  and,  from  this,  the  average 
CVK  and  VO.  The  average  separation  ZRCOH  of  the  two  pairs  of  rays  at  the 
source  end  is  computed.  If  ZRCOH  is  less  than  the  threshold  value  RCOH,  ICOH 

is  set  to  1.  Then  the  average  separation  of  the  two  pairs  of  rays  at  the  re- 

ceiver end  is  similarly  tested  and  if  it  is  found  to  be  less  than  RCOH,  ICOH 
is  increased  by  2. 

If  coherence  occurs  only  at  one  end  or  the  other,  there  will  be  two  pairs 
of  rays  to  be  processed,  and  NLOOP  is  set  to  2.  The  first  pair  selected  is  1 

and  2 if  ICOH  = 1,  and  1 and  3 if  ICOH  = 2.  The  average  value  of  VAK  and  the 

consequent  values  of  CVK  and  VO  are  computed  for  these  two  rays,  and  the  co- 
herence factor  and  resultant  coherent  intensity  are  computed  and  added  to  TRM 
as  in  the  two-ray  case. 

The  second  pair  of  rays  is  now  selected,  rays  3 and  4 if  ICOH  = 1 and  rays 
2 and  4 if  ICOH  = 2.  These  rays  are  processed  in  the  same  manner  as  the  first 
pair  and  the  resultant  intensity  is  added  to  TRM,  which  now  contains  the  com- 
bined intensity  of  all  four  rays. 

When  ICOH  = 3,  the  resultant  intensity  is  computed  in  accordance  with 
equation  (J-ll)  or  (J- 14) . These  formulas  are  identical  except  for  an  inter- 
change of  <t>^  and  (fig  (PHIA  and  PHIB)  and  for  the  method  of  calculating  g+  and 
g. . The  parameter  which  controls  the  selection  is  SN.  If  SN  = +1,  the  true 
source  is  located  at  the  ray  source  and  equation  (J-ll)  applies.  If  SN  = -1, 
the  true  source  is  located  at  the  ray  receiver  and  equation  (J-14)  applies. 

The  first  factor  in  these  equations  is  computed  and  limited,  if  necessary,  by 
the  attenuation  factor  FRACTS  in  the  same  manner  as  in  the  two-ray  case.  The 
second  factor  is  treated  in  a similar  manner  except  that  the  coefficient  in- 
volving the  g functions  is  not  present. 

In  all  the  cases  considered  the  resultant  intensity  is  stored  in  the  vari- 
able TRM.  At  the  end  of  the  frequency  loop,  TRM  is  added  to  the  contents  of 
the  appropriate  range- frequency  bin  of  the  H(K,IF)  array. 

Termination  of  Computations  Within  the  Sector 

The  criteria  for  limiting  the  number  of  ray  cycles  (NCYC)  to  be  computed 
in  a sector  are  contained  in  PLRAY  and  will  be  discussed  below.  One  of  the 
criteria  provides  for  termination  when  the  intensities  of  all  the  ray  arrivals 
in  the  sector  are  less  than  10“13,  corresponding  to  a propagation  loss  greater 
than  130  db.  The  intensity  test  is  performed  in  RAYSIIM.  Tne  maximum  intensity 
is  denoted  by  the  variable  HLIM.  Prior  to  the  beginning  of  the  big  range  loop 
IILIM  is  set  to  0.  Then  as  each  intensity  HTRM(J, JJ ,IF)  is  computed,  it  is 
tested  against  HLIM  and  if  it  is  larger,  the  value  of  IILIM  is  appropriately 
increased.  At  the  end  of  the  range  loop  HLIM  contains  the  largest  intensity 
of  all  the  ray  contributions  in  the  sector.  If  HLIM  is  less  than  10'15  and 
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the  current  sector  is  the  final  one  (IL  = NL) , the  cycle  counter  NCYC  is  set 
to  999. 

PROGRAM  PLRAY  (CONTINUED) 

Initialization  of  Ranges  and  Range  Derivatives  for  First  and  Higher  Order 
Arrivals 


After  completion  of  the  processing  of  the  zero  order  arrivals,  initial 
values  of  R(ITH,J)  and  RP(ITH,J)  are  computed  such  that  the  proper  values  for 
the  first  and  subsequent  arrival  orders  may  be  obtained  by  successively  adding 
the  full  cycle  increments  RCYC(ITH)  and  RPCYC(ITH) . Two  different  sets  of 
formulas  are  required,  depending  upon  the  value  of  MJ.  If  MJ  = 1,  only  two 
multipath  ray  types  are  involved;  if  MJ  = 2,  all  four  are  present. 

Arrival  Order  Loop;  Termination  of  Computations  in  Sector 

After  initialization,  the  program  enters  an  arrival  order  loop,  the  ar- 
rival order  number  being  indicated  by  NCYC.  At  the  beginning  of  the  loop  the 
ray  ranges  R(ITH,J)  and  range  derivatives  RP(ITH,J)  are  updated  by  adding  the 
ray  cycle  terms  RCYC(ITH)  and  RPCYC(ITH).  Before  calling  RAYSUM  to  process  the 
current  arrival  order,  a test  is  made  to  determine  whether  the  computations 
should  be  terminated. 

There  are  two  criteria  for  terminating  the  computations.  First,  in  all 
sectors  except  the  last,  no  further  computations  are  required  if  all  of  the 
ranges  R(ITH,J)  exceed  the  maximum  receiver  range  RMAX.  This  criterion  does 
not  apply  to  the  final  sector,  which  includes  bottom  bounce  propagation  at 
steep  angles  beyond  the  sector  boundary  at  30  degrees.  The  criterion  applied 
in  the  final  sector  is  based  on  intensity. 

During  the  updating  of  the  ranges,  the  variable  RM  is  set  equal  to  the 
smallest  of  all  the  R(ITH,J).  If  RM  exceeds  RMAX  and  the  current  sector  is 
not  the  final  one,  the  program  jumps  to  the  end  of  the  sector  loop,  updates 
TH1,  and  begins  to  work  on  the  next  sector.  On  the  other  hand,  if  RM  is  less 
than  RMAX  or  if  this  is  the  final  sector,  the  algebraic  signs  SIGN(J)  are 
evaluated  and  RAYSUM  is  called. 

As  mentioned  previously,  RAYSUM  computes  the  maximum  intensity  HLIM  of 
all  the  arrivals  in  the  sector,  and  if  HLIM  is  less  than  10'13  and  the  current 
sector  is  the  final  one,  NCYC  is  set  to  999.  Upon  the  return  from  RAYSUM, 

PLRAY  tests  NCYC  and  if  its  value  is  999,  the  computations  are  terminated  and 
the  program  jumps  out  of  the  sector  loop.  Otherwise  TH1  is  updated  and  the 
processing  of  the  next  sector  is  begun. 

Invocation  of  the  Surface  Duct  Model 

After  the  last  sector  has  been  processed,  the  flag  ISD  is  checked.  A 
nonzero  value  signifies  that  the  surface  duct  model  is  required  and  SUBROUTINE 
SDUCT  is  called.  A description  of  SDUCT  is  given  below. 
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Propagation  Loss 

At  the  conclusion  of  the  computations  in  RAYSUM,  RCAUST,  and  SDUCT,  the 
H(K, IF)  array  contains  the  resultant  intensities  at  all  the  ranges  and  fre- 
quencies. If  any  of  the  range- frequency  bins  is  empty,  due  to  the  absence  of 
any  rays  propagating  there,  it  is  arbitrarily  assigned  a value  HMIN,  equivalent 
to  a propagation  loss  of  999  db.  The  intensities  in  all  bins  are  then  con- 
verted to  propagation  loss  and  printed. 

Propagation  Loss  Plot 

The  plot  routine  SUBROUTINE  PLPLOT  is  called  if  IPLOT  has  been  assigned  a 
value  of  1.  This  routine  plots  the  propagation  loss  as  a function  of  range. 

SUBROUTINE  SDUCT 

This  subroutine,  together  with  the  auxiliary  routine  FUNCTION  DEPFN,  con- 
tains the  surface  duct  model  described  in  appendix  K.  The  surface  duct  model 
actually  consists  of  two  models,  a poor  duct  model  and  a good  duct  model.  The 
poor  duct  model  is  applicable  to  ducts  in  which  only  the  first  mode  of  a normal 
mode  solution  contributes  significantly  to  the  propagation.  The  good  duct 
model  applies  to  ducts  in  which  two  or  more  modes  must  be  considered. 

At  the  beginning  of  SDUCT  is  a short  section  in  which  the  beam  pattern 
(if  any)  is  computed.  The  beam  pattern  function  is  based  on  the  source  angle 
of  the  limiting  ray  to  the  bottom  of  the  duct.  Next,  the  cut-off  frequency 
FMIN  is  computed.  The  cut-off  frequency  is  a somewhat  arbitrarily  selected 
value  such  that  at  frequencies  less  than  FMIN  the  duct  is  so  leaky  that  it  can 
be  ignored.  The  remainder  of  the  subroutine  is  enclosed  within  a frequency  DO 
loop.  All  that  follows,  therefore,  applies  separately  to  each  individual  fre- 
quency. 

At  the  beginning  of  the  loop  the  frequency  (denoted  by  FR  in  SDUCT)  is 
tested  against  FMIN.  If  FR  is  less  than  FMIN,  the  program  jumps  to  the  end  of 
the  loop.  If  not,  the  next  step  is  to  compute  the  number  of  trapped  modes  EN 
(designated  by  the  symbol  n in  appendix  K) , the  leakage  and  surface  scattering 
attenuation  coefficients  ALPH1  and  ALPHS  (oi  and  a§) , and  the  parameter  VO 
(Ac0) . These  quantities  are  needed  for  both  the  poor  duct  model  and  the  good 
duct  model.  Then  EN  is  tested  against  the  limiting  value  of  1.5  to  determine 
which  of  the  two  models  is  to  be  invoked. 

Poor  Duct  Model 

The  first  step  is  the  computation  of  the  depth  functions  US  and  UR.  These 
are  computed  in  FUNCTION  DEPFN,  described  below.  Their  product  is  the  complete 
depth  function  U.  The  remainder  of  this  portion  of  the  subroutine  is  a range 
loop  (index  K)  in  which  the  intensities  are  computed  from  two  different  formu- 
las. UK  is  computed  from  the  portion  of  the  normal  mode  formula  (K- 16)  exclu- 
sive of  the  attenuation  factor,  and  HKK  is  computed  from  the  spherical  spread- 
ing formula  (K-42) , exclusive  of  the  same  attenuation  factor.  HK  and  HKK  are 
compared  and  the  larger  of  the  two  is  chosen  and  multiplied  by  the  attenuation 
and  beam  pattern  factors.  If  the  resulting  intensity  is  greater  than  the  lower 
limit  HM,  corresponding  to  130  db  propagation  loss,  it  is  added  into  the  proper 
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bin  of  the  H(K,IF)  array.  If  it  is  less  than  HM,  computations  are  terminated 
and  thw  program  jumps  to  the  end  of  the  range  loop. 

Good  Duct  Model 

This  portion  of  the  subroutine  begins  with  the  computation  of  various 
range- independent  quantities,  including  the  attenuation  coefficient  ALPH2  (02 
in  appendix  K) , the  constants  CR  and  CS  (KR  and  Kg) , and  several  other  constant 
factors  involved  in  the  intensity  formulas.  As  in  the  case  of  the  poor  duct 
model,  the  intensities  are  computed  as  a function  of  range  in  a DO  loop  over 
the  index  K.  The  intensity  HK  is  the  portion  of  Hj  exclusive  of  the  sine  fac- 
tors in  the  cylindrical  spreading  formula  (K-44)  ; similarly,  HKK  is  the  portion 
of  H2  exclusive  of  the  same  sine  factors  in  the  spherical  spreading  formula 
(K-45) . HK  and  HKK  are  compared  and  the  larger  is  selected  and  multiplied  by 
the  sine  factors.  If  there  is  a beam  pattern,  the  beam  intensity  factor,  the 
square  of  BMFAC(IF),  is  also  included.  The  resulting  intensity  is  then  checked 
against  HM  and  this  portion  of  the  routine  closes  in  the  same  manner  as  the 
poor  duct  portion. 

FUNCTION  DEPFN(Z) 

Except  for  the  difference  in  symbols,  the  coding  of  this  routine  is  a 
straightforward  application  of  the  formulas  developed  in  appendix  K.  To  aid 
in  correlating  the  FORTRAN  with  the  algebra,  several  of  the  symbols  are  listed 
in  table  I. 

SUBROUTINE  BEAM(L) 

This  subroutine  is  divided  into  two  parts.  The  first  part,  which  is  ac- 
tivated when  BEAM  is  called  from  PLRAY  with  L = 0,  calls  for  reading  the  beam 
inputs  and  printing  them.  The  second  part,  which  is  activated  when  BEAM  is 
called  from  RAYSUM,  RCAUST,  or  SDUCT  with  L = 1,  computes  the  beam  pressure 
ratios  BMFAC(IF)  corresponding  to  the  current  value  of  the  ray  tangent  VBM. 

The  first  data  card  of  part  1 contains  a set  of  integers  (NBF(IF),  IF  = 1, 
NF) . NBF(IF)  is  the  number  of  data  points  required  to  specify  the  beam  pattern 
for  the  frequency  F(IF),  the  maximum  permissible  number  being  100.  This  is 
followed  by  two  sets  of  cards  for  the  first  frequency,  then  two  sets  for  the 
second  frequency,  and  so  on.  The  first  set  of  cards  contains  the  beam  angles 
THETA(I)  in  increasing  order  from  -90  degrees  (downward)  to  +90  degrees  (up- 
ward), 10  values  per  card.  The  second  set  contains  the  corresponding  beam 
pressure  ratios  BM(I,IF)  (pressure  at  angle  0 divided  by  pressure  along  beam 
axis).  The  pressure  ratio  must  be  a real  number,  positive  if  the  pressure  is 
in  phase  with  the  pressure  along  the  axis,  and  negative  if  180  degrees  out  of 
phase.  There  is  no  provision  for  complex  ratios  corresponding  to  phases  other 
than  0 and  180  degrees.  After  the  angles  and  beam  functions  for  each  frequency 
are  read  in,  they  are  printed.  Then  the  sines  of  the  angles  are  computed  and 
stored  in  an  array  SINE (I, IF). 

The  second  portion  of  BEAM  begins  with  the  conversion  of  the  input  tangent 
function  VBM  to  a sine  function  SBM.  The  output  beam  pressure  ratios  BMFAC(IF) 
are  then  computed  in  a frequency  DO  loop  by  interpolating  in  the  SINE(I,IF)  and 
BM(I , IF)  tables. 
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TABLE  I 

FORTRAN  AND  ALCEBRAIC  SYMBOLS  IN  DEPFN 
FORTRAN  Algebraic 


B1 

all 

Cl 

a21 

C2 

a22 

C3 

a23 

D1 

a31 

D2 

a32 

F 

f 

FF 

f' 

THETA 

$ (*) 

UFAC 

A 

V 

Ac 

VO 

Aco 

VI 

AC1 

vz 

1 - 2/ 

ZA 

za 

ZAA 

za' 

Z1 

Z1 

Z2 

z2 
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COMPARISON  WITH  OTHER  MODELS 


To  compare  the  performance  of  PLRAY  with  that  of  FACT,  runs  were  made  on 
both  programs  for  a variety  of  environmental  inputs.  Runs  were  made  both  semi- 
coherently  and  incoherently.  As  a check  on  the  semicoherent  runs  of  both  pro- 
grams, comparison  runs  were  made  also  on  a locally  generated  normal  mode  pro- 
gram called  AP2  (appendix  L) , which  is  capable  of  accepting  the  same  velocity 
profile  and  bottom  loss  inputs  as  the  two  ray  programs. 

The  principal  limitation  of  a normal  mode  solution  in  deep  water  is  the 
large  number  of  modes  which  must  be  computed.  The  number  of  modes  required 
for  a complete  solution  (in  the  practical  sense)  is  approximately  equal  to  the 
number  of  half  wavelengths  in  the  ocean  depth.  For  example,  if  the  ocean  is 
15,000  feet  deep,  the  number  of  modes  required  at  100  Hz  is  approximately  600. 
Since  AP2  is  dimensioned  only  for  500  modes,  there  was  a shortage  of  modes  in 
several  of  the  runs  made  to  check  PLRAY  and  FACT.  Failure  to  compute  the 
higher  order  modes  is  equivalent  in  ray  theory  to  failure  to  compute  the 
steeper  rays.  From  the  phase  velocity  of  the  highest  order  mode  computed  it 
is  possible  to  compute  the  angle  of  the  steepest  equivalent  ray,  and  from  this 
the  minimum  range  at  which  reliable  bottom-bounce  propagation  is  computed.  In 
presenting  the  AP2  results  the  portion  of  each  curve  below  the  minimum  reli- 
able range  has  been  omitted.  It  should  be  noted,  however,  that  in  the  direct 
propagation  zone  at  very  short  ranges,  where  according  to  ray  theory  the  domi- 
nant contribution  to  the  resultant  intensity  comes  from  shallow-angle  rays 
which  propagate  directly  from  the  source  to  the  receiver,  the  absence  of  the 
higher  order  modes  does  not  introduce  appreciable  error.  For  this  reason,  in 
many  of  the  runs  which  suffer  from  a deficiency  of  modes  there  is  a short  in- 
terval (not  shown  in  the  plots)  at  the  beginning  of  the  range  scale  where  re- 
liable results  are  obtained,  the  extent  of  the  interval  depending  upon  the 
nature  of  the  velocity  profile  and  upon  the  source  and  receiver  depths.  The 
existence  of  this  interval  has  proved  helpful  in  a few  instances  in  checking 
the  performance  of  PLRAY  and  FACT. 

To  facilitate  comparison  with  the  semicoherent  predictions  of  the  two  ray 
programs  it  is  desirable  to  smooth  out  the  rapid  oscillations  of  the  fully  co- 
herent normal  mode  curves.  The  smoothing  was  accomplished  by  means  of  a 
weighted  sliding  window  containing  nine  terms^  The  losses  were  first  converted 
to  intensities.  Then  the  average  intensity  Ij  at  the  ith  range  was  computed 
from  the  formula 


I. 


i = Cri. 


.+21.  _+3I . +41 . ,+51.  +41.  .+31.  _+2I . ,+I.  J/25 

4 1-3  1-2  l-l  l l+l  1+2  1+3  1+4 


The  average  intensities  were  then  converted  back  to  db.* 


* The  actual  procedure  included  one  additional  refinement.  To  compensate  for 
the  distortion  introduced  by  inverse  square  spreading  at  short  ranges,  the 
intensities  were  multiplied  by  the  squares  of  the  respective  ranges  before 
averaging.  The  resulting  average  was  then  divided  by  the  square  of  the  cen- 
tral range  of  the  window  before  conversion  to  db. 
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The  results  for  the  semicoherent  runs  are  presented  in  the  form  of  three 
plots.  The  PLRAY  plot  is  shown  at  the  top  of  the  page,  the  FACT  plot  is  in 
the  center,  and  at  the  bottom  the  smoothed  AP2  curve  is  superimposed  on  the 
plot  of  the  original  AP2  outout. 

Both  semicoherent  and  incoherent  runs  were  made  with  eight  different  deep- 
water velocity  profiles  selected  from  a number  of  different  ocean  areas.  The 
eight  profiles  are  shown  in  figures  6 through  13.  The  bottom  loss  curves  as- 
sociated with  these  profiles  are  shown  in  figures  14  through  19.  The  curves  of 
profiles  1 and  6 (figures  14  and  18)  are  based  on  NAVAIRDEVCEN  experimental 
measurements;  the  remainder  are  FNWC  curves.  It  will  be  noted  that  the  curves 
of  figure  16  apply  to  profiles  3,  5,  and  7. 

The  runs  on  the  first  seven  profiles  were  carried  out  to  a range  of  100 
kyd.  The  runs  on  profile  8 were  limited  to  a range  of  50  kyd. 

The  propagation  loss  plots  for  most  of  the  runs  can  be  divided  into  range 
intervals  dominated  by  one  or  another  of  three  different  types  of  propagation. 
At  short  ranges  is  the  direct  propagation  zone,  dominated  by  rays  which  propa- 
gate directly  from  the  source  to  the  receiver  either  with  no  intervening  ver- 
tices or  with  a single  surface  bounce  or  upper  refractive  vertex.  These  are 
the  zero  order  arrivals.  Included  in  this  category  also  is  surface  duct  prop- 
agation out  to  the  range  at  which  the  ducted  intensity  drops  to  the  level  of 
the  bottom-reflected  intensity.  In  the  second  category  are  the  convergence 
zones,  formed  by  rays  which  are  refracted  in  the  deep  ocean  without  striking 
the  bottom.  Most  of  the  profiles  lead  to  a single  zone  within  the  range  in- 
terval of  the  computations.  In  the  third  category  is  the  bottom-bounce  region. 
To  be  precise,  this  category  should  be  broken  down  into  a series  of  regions, 
the  first  dominated  by  single  bottom-bounce  rays,  the  second  by  double  bottom- 
bounce  rays,  and  so  on.  In  most  of  the  runs  the  major  portion  of  the  plot 
consists  of  the  first  bottom-bounce  region  extending  out  to  the  convergence 
zone,  with  a small  portion  of  the  second  region  appearing  beyond  the  zone. 

PROFILE  1 

The  semicoherent  predictions  of  PLRAY  and  FACT  for  Profile  1 are  compared 
with  AP2  in  figure  20.  The  frequency  is  105  Hz  and  the  source  and  receiver 
depths  are  80  and  100  feet.  It  will  be  noted  that  there  is  remarkably  good 
agreement  among  all  three  models  at  all  ranges.  All  the  features  of  the  inter- 
ference pattern  in  the  first  bottom-bounce  region  out  to  70  kyd  are  present. 

The  agreement  between  PLRAY  and  AP2  here  and  in  the  relatively  featureless  sec- 
ond bottom-bounce  region  is  excellent.  FACT  also  agrees  very  well  except  for 
a discrepancy  in  the  ranges  at  which  the  scallops  appear  in  the  first  bottom- 
bounce  region.  The  FACT  range  scale  appears  to  be  compressed  by  about  7 per- 
cent. This  is  a phenomenon  which  has  been  observed  in  many  of  the  runs.  Some- 
times the  compression  (or  expansion)  relative  to  AP2  occurs  with  one  of  the 
ray  models,  sometimes  with  the  other,  sometimes  with  both.  The  cause  of  the 
scale  distortion  is  not  understood.  A number  of  possible  explanations  have 
been  examined  but  none  have  proved  adequate.  The  problem  will  be  discussed  in 
a later  section. 


One  curious  feature  of  figure  20  is  the  appearance  of  a convergence  zone 
at  about  76  kyd  in  both  the  PLRAY  and  FACT  curves,  whereas  there  is  not  even  a 
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FIGURE  8 - Profile  3 
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FIGURE  9 - Profile  4 
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FIGURE  11  - Profile  6 
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FIGURE  13  - Profile  8 
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FIGURE  18  - Bottom  Loss,  Profile  6 


FIGURE  19  - Bottom  Loss,  Profile  8 
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FIGURE  20  - Semicoherent  Propagation  Loss,  Profile  1;  Frequency  105  Hi, 
Source  Depth  80  Ft,  Receiver  Depth  100  Ft 
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suggestion  of  such  a phenomenon  in  the  AP2  curve.  The  presence  of  the  zone  in 
the  ray  program  outputs  is  an  example  of  a fundamental  error  of  ray  theory  and 
will  be  discussed  later. 

The  incoherent  predictions  of  PLRAY  and  FACT  are  shown  in  figure  21.  Ex- 
cept for  a minor  difference  in  the  convergence  zone,  the  two  curves  are  virtu- 
ally identical.  The  difference  in  the  convergence  zones  is  immaterial  since 
normal  mode  theory  predicts  that  no  zone  exists.  The  dip  in  both  curves  at 
about  10  kyd  is  due  to  the  sharp  rise  in  the  bottom  loss  curve  around  30  de- 
grees (figure  14).  The  compression  of  the  FACT  range  scale  relative  to  PLRAY, 
noted  previously,  is  also  visible  in  the  vicinity  of  the  dip. 

PROFILE  2 

Two  sets  of  runs  were  made  on  Profile  2,  both  at  a frequency  of  100  Hz 
and  a source  depth  of  300  feet.  The  semicoherent  results  for  the  first  set  of 
runs,  in  which  the  receiver  was  located  at  a depth  of  60  feet,  are  shown  in 
figure  22.  In  this  case  the  features  of  the  interference  pattern  in  both  bot- 
tom-bounce regions  are  reproduced  in  all  three  plots,  but  the  agreement  is  not 
as  close  as  in  the  run  on  Profile  1.  Here  the  range  scales  of  both  PLRAY  and 
FACT  are  altered  relative  to  AP2,  the  scale  of  PLRAY  being  slightly  compressed 
and  the  scale  of  FACT  expanded.  The  differences  in  range  become  very  large  in 
the  second  bottom-bounce  region. 

All  three  models  predict  a convergence  zone  with  two  peaks.  The  agreement 
on  the  second  peak  is  fairly  good,  PLRAY  being  about  2 db  too  high  and  FACT 
about  2 db  too  low.  However,  the  agreement  on  the  first  peak  is  poor,  the  nor- 
mal mode  curve  being  some  6 to  7 db  below  both  ray  curves. 

In  the  direct  propagation  zone  PLRAY  and  FACT  agree  extremely  well  with 
each  other,  but  both  deviate  from  AP2,  which  shows  a considerably  more  rapid 
increase  in  loss  in  the  interval  between  1 and  2 kyd.* 

The  incoherent  predictions  of  PLRAY  and  FACT,  plotted  in  figure  23,  are 
almost  identical  everywhere  except  in  the  convergence  zone.  It  is  curious  that 
while  there  is  little  difference  in  the  PLRAY  convergence  zone  between  the 
semicoherent  and  incoherent  runs,  there  is  considerable  difference  in  the  FACT 
runs.  The  first  peak  is  higher  and  the  dip  between  the  two  peaks  is  completely 
filled. 

The  second  set  of  semicoherent  runs,  made  for  a receiver  depth  of  200 
feet,  are  plotted  in  figure  24.  Here  the  agreement  among  the  three  models  is 
slightly  better  than  in  the  runs  at  60  feet.  In  the  bottom-bounce  regions 
PLRAY  and  FACT  yield  roughly  comparable  results  relative  to  AP2,  except  that 
the  large  dip  just  before  the  convergence  zone  is  absent  in  the  FACT  curve  and 
the  heights  of  some  of  the  peaks  do  not  agree  as  well.  FACT  shows  better 


* This  is  one  of  the  runs  for  which  the  number  of  modes  required  (approximately 
725)  is  in  excess  of  the  available  500.  The  AP2  output  is  not  plotted  in 
figure  22  at  ranges  less  than  the  minimum  reliable  range  of  12  kyd.  However, 
examination  of  the  original  computer  output  indicates  that  reliable  results 
are  obtained  in  the  direct  zone  out  to  approximately  2 kyd. 
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FIGURE  21  - Incoherent  Propagation  Loss,  Profile  1;  Frequency  105  Hz 
Source  Depth  80  Ft,  Receiver  Depth  100  Ft 
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FIGURE  22  - Semicoherent  Propagation 
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FIGURE  24  - Semi coherent  Propagation  Loss,  Profile  2;  Frequency  100  Hz 
Source  Depth  300  Ft,  Receiver  Depth  200  Ft 
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agreement  in  the  second  bottom-bounce  region  than  PLRAY.  The  range  scale  dis- 
tortions of  the  interference  patterns  are  in  the  same  direction  as  in  figure 
22,  though  somewhat  smaller  in  magnitude.  In  the  zone  itself  FACT  agrees 
closely  with  AP2;  the  PLRAY  zone  is  centered  at  too  large  a range  and  does  not 
have  the  proper  shape.  In  the  direct  zone  the  performance  of  the  three  pro- 
grams is  essentially  the  same  as  in  the  runs  at  the  60-feet  receiver  depth. 

The  incoherent  curves  for  PLRAY  and  FACT,  figure  25,  are  virtually  iden- 
tical except  in  the  convergence  zone. 

PROFILE  3 

Figure  26  shows  the  semicoherent  results  based  on  Profile  3 for  a fre- 
quency of  100  Hz,  a source  depth  of  300  feet,  and  a receiver  depth  of  60  feet. 
The  features  of  the  interference  pattern  in  the  bottom-bounce  region  can  be 
recognized  in  all  three  curves,  but  the  agreement,  except  at  ranges  below  20 
kyd,  is  not  very  good.  Here  again  the  range  scales  of  PLRAY  and  FACT  are  dis- 
torted in  opposite  directions  relative  to  AP2.  In  the  second  bottom-bounce 
region  beyond  the  convergence  zone  the  distortion  is  so  great  that  the  three 
curves  have  radically  different  appearances. 

The  agreement  in  the  convergence  zone  is  also  rather  poor.  The  AP2  zone 
exhibits  two  peaks  separated  by  a deep  valley.  PLRAY  and  FACT  show  only  a 
single  peak.  The  PLRAY  zone  has  the  proper  intensity,  but  occurs  midway  be- 
tween the  two  normal  mode  peaks.  The  FACT  zone  lies  coincident  with  the  second 
normal  mode  peak,  but  its  intensity  is  about  3 db  too  low. 

Inspection  of  the  direct  propagation  zone  reveals  a large  discrepancy  be- 
tween PLRAY  and  FACT  at  ranges  below  7 kyd.  The  FACT  curve  drops  abruptly  to 
the  bottom-bounce  level  at  a range  of  2 kyd,  whereas  the  PLRAY  curve  has  a 
gradual  slope  and  does  not  reach  the  bottom-bounce  level  before  7 kyd.  Al- 
though the  AP2  plot  of  figure  26  does  not  extend  below  12  kyd,  due  to  a defi- 
ciency of  modes,  inspection  of  the  original  output  reveals  a narrow  range  in- 
terval between  1 and  1.5  kyd  where  there  is  good  agreement  with  FACT.  A check 
of  ray  paths  shows  that  when  the  velocity  profile  is  fitted  with  curvilinear 
segments,  as  in  PLRAY,  the  family  of  surface-reflected  rays  extends  theoreti- 
cally to  infinity  and  the  region  of  significant  intensity  extends  to  about  7 
kyd,  as  may  be  seen  in  figure  26.  When  straight-line  segments  are  used,  as  in 
FACT,  the  family  of  surface-reflected  rays  is  confined  to  a finite  range  in- 
terval which  extends  to  a range  somewhere  between  5 and  5.5  kyd.  It  is  not 
clear  why  the  contributions  of  these  rays  do  not  show  up  in  the  FACT  output. 

The  incoherent  curves  for  the  same  inputs  are  shown  in  figure  27.  As  has 
been  observed  in  previous  runs,  the  agreement  in  the  bottom-bounce  regions  is 
almost  perfect.  Here  again  there  is  a large  difference  in  the  FACT  convergence 
zone  between  the  semicoherent  run  and  the  incoherent  run.  In  the  incoherent 
run  the  zone  begins  at  a shorter  range,  the  first  peak  is  higher,  and  the  sec- 
ond peak  has  disappeared.  In  PLRAY  the  zones  are  virtually  identical.  In  the 
direct  propagation  zone  the  discrepancy  noted  above  is  clearly  visible. 

A second  set  of  runs  on  Profile  3 was  made  at  a receiver  depth  of  1000 
feet,  the  semicoherent  runs  being  plotted  in  figure  28.  There  is  fair  agree- 
ment among  the  three  curves  at  short  ranges,  but  the  agreement  deteriorates 
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FIGURE  28  - Semicoherent  Propagation  Loss,  Profile  3;  Frequency  100 
ource  Depth  300  Ft,  Receiver  Depth  1000  Ft 
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with  increasing  range.  It  will  be  observed  that  the  oscillations  of  the  PLRAY 
interference  pattern  decrease  in  amplitude  and  become  smooth  at  short  ranges. 
The  FACT  oscillations,  on  the  other  hand,  retain  their  amplitude  and  choppy 
nature  and  are  in  relatively  good  agreement,  at  least  at  the  shorter  ranges, 
with  the  normal  mode  oscillations.  The  discrepancy  between  the  two  ray  models 
is  the  result  of  a difference  in  the  method  of  implementing  the  amplitude  re- 
duction technique  in  cases  of  inadequate  range  sampling  and  will  be  discussed 
in  more  detail  later.  According  to  reference  (b)  FACT  is  supposed  to  begin 
cutting  down  on  the  amplitude  when  the  number  of  range  points  per  cycle  of  the 
oscillations  falls  below  6.  At  the  shorter  ranges  of  figure  28  the  number  of 
points  per  cycle  is  less  than  6,  but  the  amplitude  is  not  being  appropriately 
reduced.  The  fact  that  in  spite  of  the  failure  to  conform  to  the  stated  algo- 
rithm, the  FACT  curve  appears  to  give  good  results  suggests  that  the  choice  of 
constants  in  the  algorithm  is  too  conservative.  It  would  probably  be  better 
to  begin  the  amplitude  reduction  at  a value  less  than  6 points  per  cycle. 

Returning  to  figure  28,  it  will  be  noted  that  neither  ray  model  reproduces 
the  normal  mode  convergence  zone  accurately.  PLRAY  agrees  well  with  the  first 
peak  but  scarcely  shows  the  second  peak  at  all.  FACT  shows  both  peaks  but 
neither  the  shapes  nor  the  heights  agree  well.  Both  programs  are  in  good 
agreement  with  AP2  on  the  trailing  edge  of  the  zone. 

In  the  direct  zone  out  to  4 kyd  all  three  programs  agree  closely.  In  the 
original  AP2  output  (not  shown  in  figure  28)  there  is  a large  dip  at  2 kyd 
which  is  more  nearly  reproduced  by  FACT  than  by  PLRAY. 

The  incoherent  curves  for  the  same  case  are  plotted  in  figure  29.  Except 
in  the  convergence  zone,  where  both  are  essentially  the  same  as  their  semico- 
herent  counterparts,  the  two  models  yield  virtually  identical  results. 

Figure  30  is  a plot  of  the  semicoherent  output  of  PLRAY  and  FACT  for  a 
frequency  of  300  Hz,  the  source  and  receiver  depths  being  300  and  1000  feet. 

AP2  was  not  run  at  this  frequency  because  of  the  excessively  large  number  of 
modes  required.  The  feature  of  primary  interest  here  is  the  behavior  in  the 
bottom-bounce  region  between  5 and  45  kyd.  It  will  be  observed  that  the  os- 
cillations of  PLRAY  start  with  essentially  zero  amplitude  and  gradually  in- 
crease to  a very  large  amplitude.  The  FACT  oscillations,  on  the  other  hand, 
retain  an  essentially  constant  amplitude  throughout.  These  curves  illustrate 
dramatically  the  difference  between  the  two  models  in  their  method  of  handling 
the  range  sampling  problem.  Both  models  employ  essentially  the  same  formula 
for  the  amplitude  reduction  as  a function  of  the  number  of  points  per  cycle. 
PLRAY  uses  the  solid  curve  of  figure  J4  (appendix  J) , and  FACT  uses  the  dashed 
cuTve.  The  difference  between  the  two  programs  is  in  the  method  of  computing 
the  number  of  range  points  per  cycle.  PLRAY  computes  the  instantaneous  value 
at  each  range  point,  based  on  an  analytical  formula.  As  may  be  seen  in  figure 
30,  the  frequency  of  the  oscillations  increases  drastically  as  the  range  de- 
creases from  45  to  5 kyd.  Strict  application  of  the  curve  of  figure  J4  results 
in  the  behavior  shown. 

FACT,  on  the  other  hand,  computes  an  average  number  of  points  per  cycle 
over  the  entire  range  interval  covered  by  the  family  of  rays  in  the  sector. 
Careful  measurement  of  the  oscillations  has  shown  that  at  ranges  below  12  kyd 
the  number  of  points  per  cycle  is  less  than  2,  so  that  if  the  algorithm  were 
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correctly  implemented,  the  amplitude  of  the  oscillations  would  be  reduced  to 
zero.  On  the  other  hand,  in  the  vicinity  of  40  kyd  the  number  of  points  per 
cycle  is  well  above  the  limiting  value  of  6,  where  no  reduction  in  amplitude 
should  be  called  for. 

It  appears  that  neither  of  the  two  curves  is  really  satisfactory.  The 
PI. RAY  curve,  even  though  it  faithfully  follows  the  rules,  looks  decidedly  un- 
reasonable. The  FACT  curve  looks  more  reasonable,  but  runs  into  sampling 
problems  at  the  shorter  ranges.  Fortunately,  very  few  runs  have  been  encoun- 
tered in  which  this  problem  has  proved  severely  troublesome.  On  the  basis  of 
the  100  Hz  case  described  above,  it  would  be  desirable  to  re-examine  the  ampli- 
tude reduction  algorithm  with  a view  to  a better  selection  of  constants.  Where 
the  oscillations  are  extremely  rapid,  as  in  the  300  Hz  case,  it  is  necessary 
to  use  a smaller  range  increment. 

The  300-Hz  run  illustrates  one  other  difference  between  PLRAY  and  FACT. 

The  difference  is  more  readily  discerned  in  the  incoherent  plots  of  figure  31. 
Comparison  of  the  curves  for  the  two  models  reveals  a noticeably  higher  loss 
for  FACT  in  the  bottom-bounce  regions  on  either  side  of  the  convergence  zone. 
This  discrepancy  is  due  to  a difference  in  the  bottom  loss  which  occurs  when 
the  internally  contained  FNWC  bottom  loss  curves  are  used.  When  the  frequency 
is  intermediate  between  the  frequencies  of  two  adjacent  curves,  PLRAY  interpo- 
lates between  the  curves  to  yield  a loss  which  varies  continuously  with  fre- 
quency. FACT  does  not  interpolate.  Instead,  it  uses  each  curve  over  a prede- 
termined interval  of  frequencies  on  either  side.  Thus,  the  loss  remains  con- 
stant with  frequency  out  to  the  interval  boundary  and  then  jumps  discontinu- 
ously  to  the  next  curve.  A frequency  of  300  Hz  lies  midway  between  the  100- 
and  500-Hz  curves.  FACT  uses  the  500  Hz  curve. 

PROFILE  4 

Three  sets  of  runs  were  made  on  Profile  4 at  a frequency  of  150  Hz  and  a 
source  depth  of  200  feet,  the  receiver  depths  being  60,  300,  and  1000  feet. 

The  semicoherent  plots  for  the  60  ft  receiver  depth  are  shown  in  figure  32. 

Over  most  of  both  bottom-bounce  regions  there  is  good  agreement  among  all  three 
models.  The  range  scale  of  the  oscillations  in  FACT  is  slightly  compressed 
relative  to  PLRAY  and  AP2  in  the  first  bottom-bounce  region,  but  the  discrep- 
ancy is  quite  small.  In  the  second  bottom-bounce  region  beyond  the  convergence 
zone  FACT  and  PLRAY  show  range  errors  of  about  1 kyd  in  opposite  directions. 

The  normal  mode  convergence  zone  exhibits  two  peaks,  the  second  being 
about  4 db  higher  than  the  first.  Only  the  first  peak  shows  up  in  the  two  ray 
programs,  and  this  peak  is  4 db  too  high.  The  gradual  tapering  off  on  the 
PLRAY  curve  on  the  far  side  of  the  zone  is  in  considerably  better  agreement 
with  AP2  than  the  FACT  curve,  which  drops  sharply. 

In  the  direct  zone  AP2  has  been  found  to  give  useful  results  out  to  a 
range  of  about  7 kyd.  The  agreement  between  PLRAY  and  AP2  in  this  region  is 
good.  FACT  also  agrees  well  out  to  4 kyd  but  begins  to  deviate  slightly  on 
the  low  loss  side  at  ranges  beyond  4 kyd.  It  should  be  noted  that  the  surface 
duct  models  of  PLRAY  and  FACT  come  into  play  in  this  run  since  the  source  and 
receiver  both  lie  within  the  275  ft  surface  duct  of  Profile  4.  The  difference 
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FIGURE  31  - Incoherent  Propagation  Loss,  Profile  3;  Frequency  300  Hz, 
Source  Depth  300  Ft,  Receiver  Depth  1000  Ft 
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in  behavior  of  the  two  ray  programs  is  probably  due  to  a difference  in  the  re- 
spective surface  duct  models. 

The  corresponding  incoherent  curves  for  this  case  are  plotted  in  figure 

33.  Except  for  the  discrepancies  noted  above,  PLRAY  and  FACT  are  in  good 
agreement. 

The  semicoherent  results  for  the  300  ft  receiver  depth  are  shown  in  figure 

34.  Here,  except  for  small  distortions  of  the  range  scales,  there  is  good 
agreement  among  the  three  programs  in  the  first  bottom-bounce  region.  The 
agreement  in  the  second  bottom-bounce  region  is  considerably  poorer,  both  ray 
programs  exhibiting  too  high  a loss  immediately  beyond  the  convergence  zone. 

The  reduction  of  the  amplitude  of  the  PLRAY  oscillations  at  ranges  below  20 
kyd  is  due  to  the  range  sampling  effect  discussed  previously.  In  the  conver- 
gence zone  itself  the  two  peaks  of  PLRAY  agree  with  the  corresponding  peaks  of 
AP2  within  about  2 db.  The  second  peak  of  FACT  occurs  at  too  large  a range 
and  between  the  two  peaks  there  is  an  excessively  deep  minimum. 

In  the  direct  propagation  zone  there  is  fairly  good  agreement  between 
PLRAY  and  AP2  (not  shown  in  figure  34).  FACT,  on  the  other  hand,  shows  a sig- 
nificantly different  trend,  dropping  to  the  bottom-bounce  level  at  about  4 kyd 
as  compared  with  7 kyd  for  PLRAY.  The  error  in  this  case  is  much  larger  and 
in  the  opposite  direction  from  that  observed  at  the  60  ft  receiver  depth.  At 
300  feet  the  receiver,  though  slightly  below  the  bottom  of  the  surface  duct, 
is  within  the  range  of  depths  covered  by  the  surface  duct  models.  Although  the 
PLRAY  and  FACT  runs  have  not  been  analyzed  to  determine  the  sources  of  the  ma- 
jor contributions  to  the  intensity,  it  is  possible  that  the  discrepancy  between 
the  two  outputs  may  be  caused  in  part  by  a difference  in  the  surface  duct  mod- 
els. 

The  incoherent  results  for  this  case  are  plotted  in  figure  35.  Here  again 
we  see  excellent  agreement  in  the  bottom-bounce  regions. 

Figure  36  shows  the  semicoherent  plots  for  the  1000  ft  receiver  depth. 

Here  the  agreement  among  the  three  curves  has  deteriorated  somewhat.  Both 
PLRAY  and  FACT  show  the  major  features  of  the  interference  pattern  reasonably 
well  out  to  a range  of  50  kyd.  The  reduction  in  amplitude  of  the  PLRAY  oscil- 
lations at  short  ranges  is  due  to  the  attenuation  introduced  by  the  range  sam- 
pling algorithm.  The  smooth  nature  of  the  oscillations  below  30  kyd  indicates 
that  only  the  single  bottom-bounce  rays  are  contributing  to  the  interference 
pattern;  the  oscillations  caused  by  the  double  bottom-bounce  rays,  which  have 
a higher  frequency,  are  completely  attenuated.  At  ranges  beyond  30  kyd  the 
double  bottom-bounce  rays  begin  to  contribute,  producing  the  minor  oscillations 
superimposed  on  the  basic  interference  pattern.  Comparison  with  AP2  shows  a 
certain  amount  of  correlation  in  the  secondary  oscillations.  The  FACT  curve 
exhibits  only  smooth  oscillations  of  uniform  amplitude  in  this  region.  Appar- 
ently the  average  number  of  points  per  cycle  computed  by  FACT  for  the  single 
bottom- bounce  rays  has  resulted  in  a partial  attenuation,  while  the  average 
number  of  points  per  cycle  for  the  double  bottom-bounce  rays  is  such  as  to 
eliminate  the  secondary  oscillations  completely. 

Beyond  50  kyd,  there  is  qualitative  agreement  between  PLRAY  and  AP2.  The 
average  levels  agree  quite  well,  but  there  is  no  one-to-one  correspondence  in 
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FIGURE  34  - Semicoherent  Propagation  Loss,  Profile  4;  Frequency  150  Hz, 
Source  Depth  200  Ft,  Receiver  Depth  300  Ft 
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the  individual  oscillations.  In  the  FACT  output  beyond  50  kyd  the  oscillations 
have  died  out  except  for  a set  of  gentle  undulations  of  very  low  amplitude  and 
long  period. 

The  normal  mode  plot  exhibits  a barely  discernible  convergence  zone  con- 
sisting of  two  peaks,  one  at  about  o4  kyd  and  the  other  at  about  72  kyd.  PLRAY 
shows  the  first  peak  at  a slightly  smaller  range  with  the  correct  amplitude 
and  the  correct  drop-off  on  the  trailing  edge,  but  fails  to  develop  the  second 
peak.  FACT  shows  both  peaks  at  ranges  which  are  2 to  3 kyd  short.  Between 
the  two  peaks  can  be  seen  the  same  strange  deep  minimum  which  occurred  in  the 
previous  run,  now  greatly  expanded  because  of  the  greater  separation  of  the 
peaks.  The  spurious  nature  of  this  feature  is  indicated  by  its  failure  to  ap- 
pear in  either  of  the  two  normal  mode  curves. 

In  the  direct  propagation  zone  there  is  again  fairly  good  agreement  be- 
tween PLRAY  and  AP2.  The  FACT  curve  drops  somewhat  too  rapidly  to  the  bottom- 
bounce  level,  though  the  discrepancy  is  smaller  than  in  the  previous  run  at 
300  feet. 

In  the  incoherent  runs,  shown  in  figure  37,  it  is  curious  that  PLRAY  ex- 
hibits a weak  second  peak  in  the  convergence  zone  at  the  same  range  as  the 
second  peak  of  FACT. 

PROFILE  5 

Runs  were  made  on  Profile  5 at  a frequency  of  100  Hz,  a source  depth  of 
300  ft,  and  two  receiver  depths,  100  and  1000  ft.  Figure  38  is  a plot  of  the 
semicoherent  results  for  the  first  receiver  depth.  There  is  reasonably  good 
agreement  among  the  three  curves  out  to  about  30  kyd,  but  rather  poor  agreement 
beyond.  The  PLRAY  convergence  zone  is  considerably  closer  to  AP2  than  that  of 
FACT.  The  curious  blip  in  the  FACT  curve  at  66  kyd  appears  to  be  the  result 
of  some  quirk  in  the  program. 

The  differences  in  the  PLRAY  and  FACT  convergence  zones  are  also  apparent 
in  the  incoherent  curves,  figure  39. 

In  the  semicoherent  runs  at  a receiver  depth  of  1000  ft,  shown  in  figure 
40,  there  is  little  similarity  in  the  curves  generated  by  the  three  programs. 

On  the  average  the  levels  predicted  by  FACT  appear  to  be  closer  to  AP2  than  the 
levels  predicted  by  PLRAY.  Reference  to  figure  10  shows  that  both  source  and 
receiver  lie  within  a sub-surface  duct  extending  to  a depth  of  about  1900  feet. 
It  appears  that  neither  ray  program  is  able  to  handle  properly  the  many  caus- 
tics which  occur  in  the  duct.  The  choppy  nature  of  the  PLRAY  curve  indicates 
that  PLRAY  had  more  difficulty  than  FACT.  The  choppy  appearance  is  also  evi- 
dent in  the  incoherent  curve  of  figure  41. 

PROFILE  6 

Figure  42  shows  the  semicoherent  run  on  Profile  6,  the  frequency  being  50 
Hz  and  the  source  and  receiver  depths  300  and  1000  feet  respectively.  Except 
for  a range  scale  compression  of  about  7 percent,  PLRAY  shows  fair  agreement 
with  AP2  out  to  the  second  peak  of  the  first  convergence  zone  at  about  45  kyd. 
Although  at  longer  ranges  the  agreement  is  rather  poor,  there  is  a good  match 
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FIGURE  39  - Incoherent  Propagation  Loss,  Profile  5;  Frequency  100  Hz, 
Source  Depth  300  Ft,  Receiver  Depth  100  Ft 
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FIGURE  41  - Incoherent  Propagation  Loss,  Profile  5;  Frequency  100  Hz, 
Source  Depth  .300  Ft,  Receiver  Depth  1000  Ft 


85 


NADC-77296-30 


at  the  first  peak  of  the  second  convergence  zone.  FACT,  on  the  other  hand, 
shows  reasonable  agreement  only  out  to  about  20  kyd.  Examination  of  the  con- 
vergence zones  reveals  very  peculiar  behavior,  with  a slow  rise  to  the  maximum 
followed  by  a precipitous  drop.  Such  behavior  suggests  that  FACT  may  have  en- 
countered a gross  error  in  applying  the  caustic  corrections  at  the  convergence 
zones.  The  same  peculiarity  is  clearly  evident  in  the  incoherent  run,  figure 
43,  indicating  that  coherence  has  nothing  to  do  with  the  problem. 

PROFILE  7 

The  runs  on  Profile  7 were  made  at  a frequency  of  100  Hz  for  a source 
depth  of  300  ft  and  receiver  depths  of  100  and  1000  ft.  The  semicoherent  prop- 
agation loss  for  the  100  ft  depth  is  plotted  in  figure  44.  Except  at  ranges 
below  20  kyd  neither  of  the  ray  programs  agrees  very  well  with  AP2.  However, 
at  these  shorter  ranges  FACT  exhibits  quite  good  agreement.  Both  curves  show 
the  same  set  of  oscillations,  although  the  FACT  scale  is  expanded  relative  to 
AP2  and  some  of  the  oscillations  have  different  peak  values.  In  the  PLRAY 
curve  the  rapid  oscillations  have  been  greatly  smoothed  at  ranges  below  10  kyd. 
The  smoothing  results  from  the  attenuation  introduced  by  the  range  sampling 
algorithm.  If  the  FACT  model  had  actually  behaved  in  accordance  with  the 
dashed  curve  of  figure  J4,  its  oscillations  would  also  have  been  greatly 
smoothed,  since  the  average  number  of  range  points  per  cycle  in  the  five  cycles 
below  10  kyd  is  only  about  3.  This  case  is  another  example  of  the  previous 
observation  (figure  30)  that  adequate  sampling  is  obtainable  when  the  number 
of  points  per  cycle  is  considerably  less  than  6. 

Both  ray  curves  agree  roughly  with  the  normal  mode  curve  in  the  vicinity 
of  the  first  convergence  zone.  However,  the  large  dip  in  the  normal  mode 
curve  at  46  kyd  does  not  appear  in  either  ray  curve.  The  second  peak  at  49 
kyd  is  present  in  PLRAY,  but  the  level  of  the  peak  is  3 db  too  high.  Beyond 
the  first  convergence  zone  neither  curve  shows  much  similarity  to  AP2 , although 
in  the  vicinity  of  the  second  zone  around  80  kyd  PLRAY  behaves  somewhat  better 
than  FACT. 

In  the  incoherent  runs,  figure  45,  PLRAY  and  FACT  agree  very  closely  ex- 
cept in  the  second  convergence  zone,  where  discrepancies  up  to  3 db  occur. 

Figure  46  shows  the  semicoherent  results  for  a receiver  depth  of  1000 
feet.  All  three  curves  are  relatively  featureless.  Except  for  the  fact  that 
the  average  levels  are  reasonably  similar,  there  is  little  agreement  among  the 
three  curves.  At  very  short  ranges  the  PLRAY  curve  drops  much  too  rapidly. 
Between  10  and  30  kyd  it  agrees  reasonably  well  with  AP2,  if  allowance  is  made 
for  the  compression  of  the  range  scale.  It  fails  to  predict  the  large  peak  at 
35  kyd  and  shows  too  high  a loss  between  33  and  45  kyd.  FACT  also  exhibits 
numerous  discrepancies,  although  on  the  whole  it  appears  to  agree  with  AP2 
somewhat  better  than  PLRAY. 

As  may  be  seen  from  the  incoherent  curves,  figure  47,  PLRAY  tends  to  pre- 
dict somewhat  higher  losses  than  FACT.  The  largest  discrepancy  occurs  at  short 
ranges . 
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FIGURE  43  - Incoherent  Propagation  Loss,  Profile  6;  Frequency  50  Hz 
Source  Depth  300  Ft,  Receiver  Depth  1000  Ft 
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FIGURE  44  - Semicoherent  Propagation  Loss,  Profile  7;  Frequency  100  Hz 
Source  Depth  300  Ft , Receiver  Depth  100  Ft 
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FIGURE  47  - Incoherent  Propagation  Loss,  Profile  7;  Frequency  100  Hz 
Source  Depth  300  Ft,  Receiver  Depth  1000  Ft 
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PROFILE  8 

Profile  8 was  selected  as  a test  case  for  the  surface  duct  models  of  PLRAY 
and  FACT,  since  it  contains  a 318- foot  surface  duct.  To  emphasize  the  effect 
of  the  duct,  a high  bottom  loss  (FNWC  Type  5)  was  arbitrarily  selected.  Fur- 
thermore, all  rays  which  penetrate  the  duct  are  bottom  limited,  so  that  there 
are  no  convergence  zones.  The  source  was  located  at  a depth  of  250  feet.  Two 
receiver  depths  were  selected,  one  at  150  feet,  near  the  center  of  the  duct, 
and  the  other  at  450  feet,  below  the  duct  but  within  the  depth  interval  covered 
by  the  PLRAY  surface  duct  model.  For  each  receiver  depth  semicoherent  runs 
were  made  at  three  frequencies,  50,  100,  and  200  Hz,  corresponding  to  duct 
depths  of  approximately  3,  6,  and  12  wavelengths  respectively. 

The  plots  at  50  Hz  for  the  shallow  receiver  depth  are  shown  in  figure  48. 
At  50  Hz  the  duct  is  so  leaky  that  the  intensity  drops  to  the  level  of  the 
bottom-bounce  propagation  within  the  first  5 kyd.  Within  this  region  it  is 
seen  that  PLRAY  agrees  closely  with  AP2,  while  the  FACT  curve  drops  somewhat 
too  rapidly.  In  the  range  interval  between  5 and  25  kyd  the  agreement  between 
PLRAY  and  AP2  is  excellent.  Beyond  25  kyd  the  two  curves  separate  widely  until 
they  merge  again  at  about  43  kyd.  The  reason  for  the  wide  divergence  is  not 
obvious,  though  it  appears  to  occur  near  the  transition  interval  between  the 
first  and  second  bottom-bounce  regions.  The  behavior  of  FACT  beyond  5 kyd, 
however,  bears  little  resemblance  to  PLRAY  and  AP2. 

The  results  at  100  Hz  are  shown  in  figure  49.  Here  the  duct  is  more  ef- 
ficient at  trapping  energy  and  the  level  of  the  bottom-bounce  propagation  is 
not  reached  until  about  10  kyd.  Except  for  a slight  distortion  of  the  range 
scale  between  20  and  30  kyd  and  an  interval  of  too  high  intensity  between  32 
and  38  kyd,  PLRAY  agrees  well  with  AP2  throughout  the  entire  plot.  FACT,  on 
the  other  hand,  shows  a radically  different  behavior.  In  the  ducted  region 
the  loss  predicted  by  FACT  is  too  high,  though  the  discrepancy  is  less  than  at 
50  Hz.  Throughout  the  bottom-bounce  region,  although  the  average  levels  are 
reasonable,  the  interference  pattern  has  an  entirely  different  character  from 
that  of  the  other  two  programs. 

At  200  Hz,  figure  50,  the  effect  of  the  duct  is  visible  out  to  a range  of 
at  least  20  kyd.  Throughout  this  range  interval  both  PLRAY  and  FACT  show 
higher  losses  than  AP2,  though  below  about  8 kyd  the  divergence  of  PLRAY  is 
somewhat  smaller.  The  behavior  of  the  AP2  curve  at  ranges  below  5 kyd  appears 
to  be  questionable,  due  to  a shortage  of  modes.  From  15  to  30  kyd  PLRAY  shows 
the  same  major  features  as  AP2,  though  the  losses  tend  to  be  3 db  too  high. 
Beyond  30  kyd  there  is  a wide  divergence,  with  the  PLRAY  curve  continuing  to 
drop  in  large  irregular  oscillations  and  the  AP2  curve  leveling  out  at  values 
between  85  and  90  db.  FACT  shows  none  of  the  features  of  the  other  two  pro- 
grams, though  the  average  level  appears  to  follow  the  same  trend  as  PLRAY. 

The  50-Hz  curves  for  the  450  ft  receiver  depth  are  plotted  in  figure  51. 

At  very  short  ranges  below  about  3 kyd  both  ray  programs  agree  fairly  well  with 
AP2.  From  3 to  8 kyd  PLRAY  and  FACT  agree  with  each  other,  but  both  curves 
depart  widely  from  AP2,  which  drops  much  more  steeply  with  increasing  range. 

In  the  bottom-bounce  region  out  as  far  as  28  kyd  PLRAY  closely  follows  the 
normal  mode  curve.  The  two  curves  diverge  widely  in  the  interval  between  28 
and  45  kyd,  but  come  together  again  between  45  and  50  kyd.  The  behavior  of 
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FIGURE  51  - Semicoherent  Propagation  Loss,  Profile  8;  Frequency  50  Hz 
Source  Depth  250  Ft,  Receiver  Depth  450  Ft 
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FACT  beyond  7.5  kyd  is  exceedingly  strange.  It  shows  none  of  the  phase  oscil- 
lations characteristic  of  the  semicoherent  mode  of  ray  summation. 

At  100  Hz,  figure  52,  FACT  follows  AP2  quite  closely  out  to  the  break  in 
the  slope  at  7.5  kyd.  PLRAY  shows  a similar  trend,  but  the  propagation  loss 
is  about  3 db  less.  In  the  interval  from  10  to  25  kyd  there  is  fair  agreement 
between  PLRAY  and  AP2,  though  not  as  good  as  at  50  Hz.  Beyond  25  kyd  the  two 
curves  depart  widely  except  for  a short  interval  between  45  and  50  kyd.  As 
before,  FACT  shows  no  interference  pattern  at  all.  Even  more  curiously,  the 
100  Hz  FACT  curve  is  virtually  identical  to  the  50  Hz  curve.  For  some  reason 
the  semicoherent  predictions  of  FACT  fail  to  show  phase  effects  and  appear  to 
be  independent  of  frequency  between  50  and  100  Hz. 

The  200  Hz  curves  are  plotted  in  figure  53.  At  short  ranges  the  FACT 
predictions  are  virtually  the  same  as  at  50  and  100  Hz.  PLRAY  shows  the  same 
general  trend,  but  fluctuates  about  FACT  with  a maximum  difference  of  about  3 
db.  Beyond  10  kyd  the  behavior  of  PLRAY  relative  to  AP2  is  similar  to  the  be- 
havior at  lower  frequencies,  except  that  the  interference  patterns  no  longer 
agree  in  detail.  The  agreement  deteriorates  progressively  with  increasing  fre- 
quency. The  FACT  curve  beyond  the  kink  at  7.5  kyd  has  the  same  appearance  as 
at  the  lower  frequencies,  except  for  a trend  toward  increasingly  greater  loss 
with  increasing  range.  The  loss  at  50  kyd  is  about  7 db  higher.  The  absence 
of  an  interference  pattern  at  all  three  frequencies  and  the  identity  of  the  re- 
sults at  50  and  100  Hz  suggest  the  existence  of  a flaw  in  the  FACT  program. 

CUSPED  CAUSTICS 

To  investigate  the  behavior  of  PLRAY  in  the  presence  of  a cusped  caustic, 
a semicoherent  run  was  made  at  100  Hz  on  Profile  3 with  the  source  and  receiver 
both  at  300  feet.  The  results  are  shown  in  the  top  graph  of  figure  54. 

Clearly  the  use  of  the  conventional  smooth  caustic  correction  procedure  in 
this  case  leads  to  a gross  error,  the  predicted  loss  at  the  peak  of  the  con- 
vergence zone  being  only  54  db.  The  corresponding  curve  for  FACT  is  plotted 
on  the  bottom  graph  of  figure  54,  where  the  loss  at  the  convergence  zone  peak 
is  72  db. 

As  a check  on  the  sensitivity  of  the  error  to  the  difference  in  depth  be- 
tween source  and  receiver,  runs  were  made  on  PLRAY  with  the  source  at  300  feet 
and  the  receiver  at  various  depths  between  280  and  320  feet.  The  heights  of 
the  convergence  zone  peaks  were  measured  and  are  plotted  as  a function  of  re- 
ceiver depth  in  the  solid  curve  of  figure  55.  It  is  seen  that  there  is  only 

a very  gradual  rise  in  the  peak  as  the  limiting  depth  of  300  feet  is  approached 

from  either  direction,  until  the  difference  in  depth  is  reduced  to  one  foot, 
the  smallest  difference  for  which  runs  were  made.  It  is  thus  seen  that  the 
region  of  serious  error  is  very  narrow. 

Companion  runs  were  made  on  FACT  and  AP2  at  receiver  depths  of  280,  290, 
295,  300,  305,  310,  and  320  feet.  The  peak  values  of  the  FACT  convergence 
zone  are  plotted  in  the  dashed  curve  of  figure  55.  The  corresponding  peak 
values  for  AP2,  after  smoothing  by  a 5-point  weighted  sliding  window,  are  plot- 
ted in  the  dot-dash  curve.  The  selection  of  5 points  for  the  window  is  a com- 

promise between  the  unsmoothed  output,  which  contains  a local  spike  at  the 
maximum  point,  and  the  9-point  smoothing  which  introduces  too  much  rounding  of 
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FIGURE  52  - Semicoherent  Propagation  Loss,  Profile  8;  Frequency  100  Hz, 
Source  Depth  250  Ft,  Receiver  Depth  450  Ft 
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FIGURE  54  - Semicoherent  Propagation  Loss,  Profile  3,  Showing  Effect 
of  Cusped  Caustic;  Frequency  100  Hz 
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FIGURE  55  - Height  of  Convergence  Zone  Peak  as  a Function  of  Receiver  Depth 
Profile  3;  Frequency  100  Hz 


NADC-77296-30 


the  narrow  convergence  zone  peak.  A plot  of  the  unsmoothed  data  would  raise 
the  AP2  curve  of  figure  55  by  about  1 db. 

Inspection  of  figure  55  reveals  that  the  normal  mode  curve  is  essentially 
flat  over  the  entire  plot,  the  maximum  variation  being  only  0.5  db.  At  depths 
less  than  290  feet  and  greater  than  310  feet  FACT  exhibits  the  same  tendency 
as  PLRAY  with  even  higher  zone  peaks.  However,  between  290  and  310  feet  the 
curve  drops  below  the  AP2  curve,  suggesting  that  the  cusped  caustic  correction 
of  FACT  produces  a slight  over-correction. 

To  avoid  the  gross  error  resulting  from  the  lack  of  a cusped  caustic  cor- 
rection, PLRAY  has  been  modified  to  insert  a provision  to  move  the  source  up- 
ward and  the  receiver  downward,  when  necessary,  to  insure  a depth  separation 
of  at  least  2 feet.  When  this  is  done,  a diagnostic  is  printed.  The  effect 
of  the  modification  may  be  seen  in  the  center  plot  of  figure  54,  in  which  the 
source  is  at  299  feet  and  the  receiver  at  301  feet.  The  2-foot  depth  separa- 
ration  has  reduced  the  convergence  zone  to  reasonable  proportions,  within  less 
than  3 db  of  AP2,  without  altering  the  curve  perceptibly  at  ranges  well  away 
from  the  cusped  caustic. 


DISCUSSION 

PLRAY  has  been  compared  with  FACT  in  21  different  cases  involving  8 dif- 
ferent velocity  profiles.  In  most  of  the  cases  comparisons  have  been  made  for 
both  the  semicoherent  and  the  incoherent  types  of  ray  summation.  In  most  of 
the  semicoherent  cases  the  two  ray  program  outputs  have  been  checked  against 
the  output  of  the  AP2  normal  mode  program  which  provides  a more  accurate  rep- 
resentation of  the  sound  field  because  it  is  an  exact  solution  of  the  wave 
equation.  In  many  of  the  runs  the  normal  mode  program  was  not  capable  of  pro- 
ducing useful  results  at  short  ranges  in  the  bottom-bounce  region  where  the 
number  of  modes  required  was  in  excess  of  the  maximum  available  number  of  500. 
The  portions  of  the  AP2  propagation  loss  plots  below  the  minimum  reliable  range 
have  been  omitted.  However,  in  most  of  these  cases  AP2  produced  useful  results 
in  the  direct  propagation  zone  at  very  low  ranges  and  although  not  shown  on  the 
plots,  these  results  have  provided  a check  on  the  performance  of  PLRAY  and 
FACT. 

INCOHERENT  RUNS 

It  has  been  observed  in  the  incoherent  runs  that  when  using  the  same  bot- 
tom loss  inputs  PLRAY  and  FACT  consistently  show  excellent  agreement  in  the 
bottom-bounce  regions.  The  one  case  (figure  31)  in  which  a significant  dif- 
ference was  observed  was  a run  at  a frequency  of  300  Hz.  The  discrepancy  is 
due  to  the  manner  in  which  the  FNWC  bottom  loss  curves  are  implemented.  In  the 
low  frequency  region  FNWC  curves  are  provided  at  frequencies  of  100,  500,  and 
1000  Hz.  At  all  frequencies  below  100  Hz  the  100  Hz  curve  is  used  in  both 
programs.  Between  100  and  500  Hz  and  between  500  and  1000  Hz  PLRAY  interpo- 
lates linearly  between  the  appropriate  pair  of  curves.  FACT,  on  the  other 
hand,  applies  the  100  Hz  curve  to  all  frequencies  up  to  150  Hz,  the  500  Hz 
curve  to  all  frequencies  between  150  Hz  and  700  Hz,  and  the  1000  Hz  curve  to 
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all  frequencies  between  700  and  1000  Hz;  thus,  there  are  discontinuous  jumps 
in  the  bottom  loss  at  150  and  700  Hz. 

The  behavior  of  PLRAY  and  FACT  in  the  direct  propagation  zone  (the  region 
covered  by  the  direct  and  surface-reflected  shallow  rays)  is  much  the  same  in 
both  the  incoherent  and  semicoherent  runs  and  will  be  discussed  below.  With  a 
few  exceptions  the  same  is  true  of  the  convergence  zones. 

SEMICOHERENT  RUNS  - BOTTOM- BOUNCE  INTERFERENCE  PATTERNS 

In  the  semicoherent  runs  the  phase  interference  patterns  are  manifested 
chiefly  in  the  bottom-bounce  regions.  Since  the  semi coherence  technique  con- 
siders only  the  coherence  between  the  rays  of  each  individual  arrival  order 
and  treats  the  resultant  intensities  of  the  various  arrival  orders  incoherent- 
ly, the  interference  patterns  tend  to  be  smoothly  scalloped.  The  normal  mode 
output  is  fully  coherent  and  contains  a fine  structure  which  results  from  in- 
terference between  different  arrival  orders.  Comparison  with  the  PLRAY  and 
FACT  outputs  is  facilitated  by  the  smoothing  technique  which  filters  out  most 
of  the  fine  structure.  It  will  be  noted,  however,  that  the  smoothing  proce- 
dure produces  excessive  rounding  of  the  sharp  dips  at  the  points  of  destruc- 
tive interference  between  the  scallops.  At  such  points  the  original  unsmoothed 
AP2  curves  provide  a more  realistic  comparison  with  PLRAY  and  FACT. 

In  most  of  the  runs  the  individual  features  of  the  AP2  smoothed  interfer- 
ence pattern  are  clearly  recognizable  in  both  the  PLRAY  and  FACT  curves.  In 
several  of  the  runs  the  agreement  is  excellent.  The  agreement  tends  in  general 
to  deteriorate  with  increasing  source  and/or  receiver  depth,  with  increasing 
range,  and  with  increasing  frequency.  The  deterioration  with  increasing  depth 
is  undoubtedly  due  to  the  approximate  nature  of  the  coherence  computation, 
which  is  based  on  the  assumption  that  between  the  depth  of  interest  (source  or 
receiver)  and  the  surface  the  direct  and  surface-reflected  rays  are  parallel 
straight  lines.  As  the  source  or  receiver  becomes  deeper,  the  rays  depart  more 
from  parallelism  and  more  distortion  is  introduced  by  ray  refraction.  Similar 
considerations  apply  to  the  effect  of  increasing  range.  At  short  ranges  where 
the  angles  are  large,  the  rays  are  approximately  straight  lines.  As  the  range 
increases,  the  rays  become  more  nearly  horizontal  and  refraction  effects  be- 
come important.  These  effects  are  aggravated  by  an  increase  in  the  frequency, 
which  shortens  the  wavelength  and  increases  the  phase  differences. 

In  comparing  the  interference  patterns  of  the  three  programs  it  has  been 
observed  in  a number  of  instances  that  although  the  patterns  agree  quite  well 
in  the  levels  and  relative  locations  of  the  individual  features,  the  absolute 
ranges  do  not  agree.  The  effect  appears  as  an  expansion  or  contraction  of  the 
range  scale.  The  most  likely  source  of  the  range  distortion  was  at  first 
thought  to  be  the  difference  in  the  methods  employed  in  the  three  programs  in 
fitting  the  velocity  profile.  FACT  uses  linear  segments;  PLRAY  uses  curvilin- 
ear segments  which  join  with  continuous  slope;  AP2  uses  segments  in  which  the 
square  of  the  index  of  refraction  (or  l/c^)  varies  linearly  with  depth.  As  a 
check  on  this  hypothesis  runs  on  Profile  1 covering  a portion  of  the  bottom- 
bounce  region  were  made  with  the  old  NADC  Ray-Tracing  Program  (reference  (a)) 
using  both  linear  and  curvilinear  segments,  since  that  program  offers  the  op- 
tion of  either  type  of  curve  fit.  The  results  were  negative.  The  change  from 
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curvilinear  to  straight-line  segments  produced  virtually  no  change  in  the  in- 
terference pattern. 

A second  possible  source  of  the  range  distortion  is  the  difference  between 
PLRAY  and  FACT  in  the  method  of  computing  the  coherent  intensity.  Where  the 
four  bottom- bounce  rays  of  a given  arrival  order  are  involved,  PLRAY  uses  the 
average  range  of  all  four  rays,  whereas  FACT  uses  only  the  first  (down-up)  ray. 
However,  a check  on  the  range  difference  in  the  run  with  Profile  1 revealed 
that  the  effect  was  too  small  in  magnitude.  Moreover,  this  effect  would  always 
be  in  the  same  direction,  whereas  errors  in  both  directions  have  been  observed 
in  the  various  runs.  In  addition,  this  explanation  does  not  account  for  the 
discrepancies  between  PLRAY  and  AP2.  A third  possible  source  of  error  is  the 
approximation  made  by  FACT  in  replacing  the  actual  range-angle  curves  with 
parabolas  fitted  to  three  points.  It  has  been  found  that  this  approximation 
is  capable  of  introducing  both  range  and  intensity  errors.  However,  it  can 
have  no  bearing  on  the  behavior  of  PLRAY  relative  to  AP2. 

The  source  of  the  range  distortion  has  not  yet  been  determined.  It  is 
clearly  profile-dependent  and,  relative  to  AP2,  it  appears  to  affect  PLRAY  and 
FACT  in  opposite  directions.  With  profiles  having  appreciable  surface  ducts 
(Profiles  1,  4,  and  8)  there  seems  to  be  a tendency  for  PLRAY  to  show  an  ex- 
pansion and  FACT  a contraction,  while  with  the  other  profiles  there  seems  to 
be  a tendency  in  the  other  direction.  Finally,  it  should  be  noted  that  the 
range  scale  distortion  affects  only  the  bottom-bounce  region;  it  appears  to 
have  no  effect  (or  at  least  no  correlated  effect)  on  the  location  of  the  con- 
vergence zones. 

If  allowance  is  made  for  the  distortion  of  the  range  scales,  there  is  fair 
to  good--and  occasionally  excellent--agreement  in  the  interference  patterns  be- 
tween AP2  and  the  two  ray  programs,  out  to  at  least  moderate  ranges,  in  most 
of  the  runs.  On  the  whole,  neither  of  the  programs  performed  significantly 
better  than  the  other.  There  are,  however,  two  notable  exceptions.  First,  in 
the  runs  on  Profile  8,  with  the  source  in  the  surface  duct,  FACT  shows  anoma- 
lous behavior  at  all  three  frequencies.  In  the  first  set  of  runs  with  the  re- 
ceiver in  the  duct  the  interference  pattern  generated  by  FACT  differs  radically 
from  those  of  PLRAY  and  AP2.  Even  more  anomalously,  in  the  second  set  of  runs, 
with  the  receiver  below  the  duct,  FACT  exhibits  no  interference  pattern  at  all. 
Moreover,  FACT  generated  the  same  output  at  100  Hz  as  at  50  Hz. 

Secondly,  in  the  two  runs  (figure  40  and  46)  in  which  both  source  and  re- 
ceiver were  located  near  the  axis  of  a major  sound  channel,  both  PLRAY  and 
FACT  failed  to  show  reasonable  agreement  with  AP2.  The  first  of  these  runs 
was  based  on  a profile  from  the  Iceland  area.  Profile  5 contains  a very  deep 
sub-surface  channel  extending  down  to  1600  feet,  with  its  axis  at  485  feet. 

The  source  was  located  at  300  feet,  well  into  the  channel,  and  the  receiver 
was  located  below  the  channel  axis  at  1000  feet.  In  such  a situation  the  dom- 
inant contribution  to  the  propagation  tends  to  come  from  the  channel  and  con- 
tains many  more  caustics  than  are  normally  encountered  when  no  such  channel  is 
present.  It  appears  that  neither  of  the  programs  did  an  adequate  job  of  han- 
dling the  caustics.  However,  as  may  be  seen  in  figures  40  and  41,  FACT  per- 
formed considerably  better  than  PLRAY.  The  second  case  of  poor  agreement  in- 
volves Profile  7,  which  was  taken  from  the  Mediterranean  Sea.  The  region  of 
minimum  sound  speed  in  Profile  7 is  very  broad  and  flat,  ranging  from  about 
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250  feet  to  1600  feet  and  encompassing  both  the  source  at  300  feet  and  the  re- 
ceiver at  1000  feet.  In  addition  to  the  caustic  problem,  which  does  not  appear 
to  have  been  so  serious  in  this  case,  there  is  another  factor  which  adversely 
affects  the  coherence  computations.  The  sound  speed  at  the  surface  is  more 
than  80  ft/sec  higher  than  at  either  the  source  or  receiver  and  in  the  inter- 
vening region  there  is  a very  strong  negative  gradient.  The  resulting  bending 
of  the  rays  may  be  expected  to  introduce  significant  errors  into  the  assumption 
of  straight-line  propagation  on  which  the  coherence  computations  are  based. 

The  range  sampling  problem  has  been  discussed  in  detail  in  connection  with 
the  runs  on  Profile  3 and  will  only  be  summarized  here  for  the  sake  of  com-- 
pleteness.  In  FACT,  whenever  the  number  of  range  points  per  cycle  of  the  in- 
terference oscillations  is  too  small  to  provide  adequate  sampling,  the  ampli- 
tude of  the  oscillations  is  reduced  by  application  of  an  attenuation  factor. 

The  attenuation  factor  varies  from  the  full  value  of  1 at  6 points  per  cycle 
to  0 at  2-2/3  points  per  cycle.  PLRAY  uses  a smooth  curve  which  closely  ap- 
proximates the  linear  segments  of  FACT.  The  two  programs  differ,  however,  in 
the  method  of  computing  the  number  of  points  per  cycle.  PLRAY  computes  in- 
stantaneous values  from  an  analytical  formula,  whereas  FACT  computes  an  aver- 
age value  based  on  the  total  range  interval  covered  by  all  the  rays  of  the 
sector.  In  cases  where  the  number  of  points  per  cycle  varies  drastically  from 
one  end  of  the  range  interval  to  the  other,  the  attenuation  factors  applied  by 
FACT  deviate  drastically  from  the  factors  applied  by  PLRAY. 

Comparison  of  the  FACT  output  with  AP2  at  short  ranges  where  the  number 
of  points  per  cycle  is  far  less  than  6 and  yet  the  amplitudes  of  the  oscilla- 
tions are  not  appropriately  attenuated,  has  led  to  the  conclusion  that  the 
values  of  6 and  2-2/3  are  too  high  and  should  be  reduced.  Further  investiga- 
tion is  needed  to  determine  the  proper  values. 

CONVERGENCE  ZONES 

In  general  there  is  quite  good  agreement  among  the  three  programs  on  the 
locations  of  the  convergence  zones  but  much  poorer  agreement  on  the  shapes  of 
the  zones  and  the  heights  of  the  peaks.  In  one  or  two  runs  (see  figures  34 
and  38)  PLRAY  gave  better  performance  in  relation  to  AP2  than  FACT;  in  others 
(see  figure  24)  FACT  was  definitely  superior.  Occasionally  (see  figure  44) 
both  ray  programs  performed  reasonably  well.  However,  in  the  majority  of  the 
runs  neither  PLRAY  nor  FACT  agreed  well  enough  with  AP2  to  make  it  meaningful 
to  attempt  to  rate  one  of  them  as  being  better  than  the  other. 

There  are  a few  runs  which  merit  special  comment.  In  the  run  on  Profile 
1 (figure  20)  the  convergence  zone  shown  by  both  PLRAY  and  FACT  is  not  present 
in  the  AP2  output.  The  appearance  of  the  zones  indicates  a basic  error  of  ray 
theory.  Reference  to  figure  6 shows  that  the  depth  excess  in  Profile  1 is  less 
than  100  feet.  This  means  that  the  energy  which  forms  the  convergence  zone  is 
confined  to  an  extremely  narrow  bundle  of  rays.  Whereas  ray  theory  assumes 
that  the  bundle  remains  intact  at  all  ranges,  the  normal  mode  output  shows  that 
there  is  actually  sufficient  leakage  by  diffraction  to  dissipate  the  energy 
before  it  reaches  the  convergence  zone  range.  In  the  runs  on  Profile  4 it  will 
be  noted  that  as  the  receiver  depth  is  increased  (figures  34  and  36)  FACT  de- 
velops a very  deep  valley  between  the  two  convergence  zone  peaks.  The  valley 
does  not  appear  in  either  PLRAY  or  AP2.  In  the  runs  on  Profile  6,  figures  42 
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and  43,  the  FACT  convergence  zones  have  a very  unrealistic  appearance.  The 
fact  that  there  is  fair  agreement  between  PLRAY  and  AP2  in  the  zone  suggests 
that  FACT  has  run  into  a problem,  probably  involving  the  caustic  correction. 

In  the  runs  on  Profiles  2,  3,  and  5,  figures  22  and  23,  26  and  27,  and  38 
and  39,  there  is  a noticeable  difference  in  the  FACT  convergence  zones  between 
the  semicoherent  and  incoherent  cases,  but  little  or  no  difference  in  the  PLRAY 
zones.  In  the  incoherent  run  of  figure  27  the  leading  edge  of  the  FACT  zone 
occurs  about  one  kiloyard  sooner  and  the  first  peak  is  about  2 db  higher  than 
in  the  semicoherent  run.  The  second  peak  has  disappeared  altogether.  In  the 
case  of  Profile  5 the  intensity  throughout  the  bulk  of  the  semicoherent  zone, 
figure  38,  is  3 to  5 db  higher  than  the  intensity  in  the  incoherent  zone,  fig- 
ure 39. 

DIRECT  ZONES 

The  primary  interest  in  the  direct  zone  propagation  is  focused  on  those 
profiles  (1,  4,  and  8)  which  contain  surface  ducts  of  appreciable  thickness. 

In  the  other  profiles,  with  one  exception,  there  is  reasonably  good  agreement 
among  the  three  programs.  The  exception  is  the  run  of  figure  26  and  has  been 
discussed  earlier. 

The  surface  duct  of  Profile  1 is  rather  poor  at  the  frequency  of  105  Hz, 
the  number  of  trapped  modes  (see  equation  K-13  of  appendix  K)  being  only  0.5. 

As  a result,  the  intensity  drops  very  rapidly  in  the  first  few  kiloyards.  As 
may  be  seen  in  figure  20,  PLRAY  and  FACT  show  essentially  the  same  rate  of 
fall-off.  Beyond  3 kyd,  where  the  normal  mode  curve  begins,  there  is  good 
agreement  also  with  AP2. 

The  duct  of  Profile  4 at  150  Hz  is  slightly  better,  the  trapped  mode  pa- 
rameter being  0.7.  The  direct  propagation  zone  extends  to  about  7 kyd.  As 
pointed  out  previously,  there  is  fair  to  good  agreement  between  PLRAY  and  AP2 
throughout  this  region  in  the  runs  at  all  receiver  depths,  figures  32,  34,  and 
36.  In  the  first  run  FACT  shows  a tendency  to  deviate  slightly  in  the  direc- 
tion of  smaller  loss.  In  the  second  and  third  runs  FACT  deviates  in  the  direc- 
tion of  higher  loss,  rather  badly  so  in  the  second  run. 

The  trapped  mode  parameter  of  Profile  8 has  values  of  0.44,  0.63,  and 
1.05  at  frequencies  of  50,  100,  and  200  Hz.  Although  all  three  values  lie  be- 
low the  threshold  at  which  the  good  duct  model  of  PLRAY  applies,  a value  in 
excess  of  1 results  in  faiTly  good  duct  propagation,  as  is  evident  in  figure 
50.  In  the  runs  in  which  the  receiver  was  located  within  the  duct,  figures 
48  through  50,  PLRAY  was  found  to  agree  well  with  AP2  at  the  two  lower  frequen- 
cies while  FACT  predicted  a loss  that  was  somewhat  too  large.  At  200  Hz,  both 
PLRAY  and  FACT  predicted  too  high  a loss  relative  to  AP2 , FACT  being  a little 
closer  than  PLRAY.  The  results  of  the  runs  for  the  case  in  which  the  receiver 
was  below  the  duct  were  quite  variable.  At  50  Hz,  figure  51,  PLRAY  and  FACT 
agree  with  each  other  but  deviate  from  AP2,  which  shows  an  appreciably  higher 
loss  in  the  direct  zone.  At  100  Hz,  figure  52,  FACT  agrees  with  AP2  while 
PLRAY  shows  too  low  a loss.  At  200  Hz,  figure  23,  no  two  of  the  curves  agree 
very  well,  though  all  three  show  essentially  the  same  trend. 
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It  is  clear  that  the  surface  duct  models  of  both  PLRAY  and  FACT  are  a ma- 
jor improvement  over  ray  theory  in  predicting  acoustic  propagation  in  a sur- 
face duct.  Both  models  predict  approximately  the  correct  leakage  of  energy 
out  of  the  duct  as  a function  of  frequency.  With  the  source  and  receiver  both 
in  the  duct,  PLRAY  appears  to  give  results  which  agree  somewhat  more  closely 
with  normal  mode  predictions  than  FACT.  When  the  source  or  receiver  is  slight- 
ly below  the  bottom  of  the  duct,  the  results  appear  to  be  more  variable  with 
neither  model  agreeing  particularly  well  with  the  normal  mode  results.  In  this 
case,  of  course,  the  only  contribution  from  the  duct  comes  from  the  leakage  of 
energy  out  of  the  duct.  Considering  the  slopes  of  the  propagation  loss  curves 
as  a function  of  frequency  it  is  seen  that  the  frequency  variation  has  far  less 
effect  when  the  source  or  receiver  is  below  the  duct  than  when  both  source  and 
receiver  are  in  the  duct. 

In  the  AP2  curves  there  is  a noticeable  change  in  the  slope  between  50 
and  100  Hz  but  very  little  change  between  100  and  200  Hz.  In  the  PLRAY  curves 
no  consistent  trend  can  be  discerned,  since  the  50  and  200  Hz  curves  are  es- 
sentially the  same  out  to  10  kyd,  whereas  the  100  Hz  curve  shows  a considerably 
gentler  slope.  The  FACT  curves  are  suspect  because  of  the  identity  of  the 
curves  throughout  their  entire  length  at  50  and  100  Hz  and  the  absence  of  an 
interference  pattern  in  all  three  curves. 

CUSPED  CAUSTICS 

Lacking  a cusped  caustic  correction,  PLRAY  gives  clearly  unacceptable  re- 
sults when  the  source  and  receiver  are  at  the  same  depth,  as  is  evident  in 
figure  54.  The  error  is  manifested  as  an  erroneous  caustic  correction  which 
results  in  an  excessively  wide  convergence  zone  with  an  erroneously  high  peak 
intensity.  However,  the  error  is  extremely  sensitive  to  the  difference  in 
depth  between  the  source  and  receiver.  In  the  sample  runs  on  Profile  3 a sep- 
aration of  only  two  feet  reduced  the  level  of  the  zone  peak  by  15  db  to  within 
less  than  3 db  of  the  values  predicted  by  the  normal  mode  program. 
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APPENDIX  A 
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BLANK  COMMON 


Symbol 

Defined 

Definition 

ALPHA (IF) 

PLRAY 

Attenuation  coefficient  of  water  (dB/kyd) 

B LOSS ( I TH , I F ) 

I NCR 

Bottom  reflection  loss  (dB) 

BMFAC(IF) 

BEAM 

Beam  loss  factor  (pressure  ratio  relative  to  beam  axis) 

BOTL(IF) 

BOTTOM 

Bottom  reflection  loss  (dB) 

CA(IA) 

INSERT 

Sound  speed,  augmented  velocity  profile  (ft/sec) 

CAA 

INSERT 

Sound  speed  at  ray  source  (ft/sec) 

CB(IB) 

PLRAY 

Sound  speed,  original  velocity  profile  (ft/sec  or  m/sec) 

CBB 

INSERT 

Sound  speed  at  ray  receiver  (ft/sec) 

CL(IL) 

SECLIM 

Sound  speed  at  sector  boundary  (ft/sec) 

CSP(IB) 

CURFIT 

Parameter  of  V.P.  profile  segment 

CV 

PLRAY 

Ray  vertex  velocity  (ft/sec) 

DR 

PLRAY 

Range  increment  (yd) 

F (IF) 

PLRAY 

Frequency  (Hz) 

F 1 3 ( I F ) 

PLRAY 

Cube  root  of  frequency 

GP(IB) 

CURFIT 

Parameter  of  V.P.  profile  segment 

H(K,IF) 

RAY SUM 
RCAUST 
SDUCT 

Resu’tant  multipath  intensity  (converted  to  propagation 
loss  in  PLRaY) 

HC(KC, J , IF) 

RCAUST 

Temporary  intensity  array 

IAA 

INSERT 

Index  of  ray  source  depth,  augmented  V.P. 

I AD 

INSERT 

Index  of  bottom  of  surface  duct,  augmented  V.P. 

IBB 

INSERT 

Index  of  ray  receiver  depth,  augmented  V.P. 

I BEAM 

PLRAY 

Beam  control  parameter 

IBOT 

PLRAY 

Flag  indicating  whether  ray  hits  bottom 

IBP 

PLRAY 

Bottom  type  parameter 

ICAUST(J) 

RAYSUM 

Caustic  type  (min.  range,  max.  range,  none) 

IDBT(4) 

PLRAY 

Title  of  Run 

IDF (I A) 

INSERT 

Identification  function  relating  IA  to  IB 

I DUCT 

PLRAY 

Index  of  bottom  of  surface  duct  (original  V.P.) 

IFMAX 

PLRAY 

Index  of  largest  frequency 

I LB 

SECLIM 

Value  of  IL  in  first  sector  in  which  rays  hit  bottom 

IliS 

SECLIM 

Value  of  IL  in  first  sector  in  which  rays  hit  surface 
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I PR 

PLRAY 

Diagnostic  print  control  parameter 

I PROF 

PLRAY 

Environmental  input  control  parameter 

ISD 

INSERT 

Flag  determining  implementation  of  surface  duct  model 

ISP 

PLRAY 

Surface  parameter  (sea  state) 

I SURF 

PLRAY 

Flag  indicating  whether  ray  hits  surface 

KA(J.IF) 

RCAUST 

Equivalent  to  K2  (RCAUST) 

KB(J, IF) 

RCAUST 

Equivalent  to  K3  (RCAUST) 

KMAX 

PLRAY 

Index  of  maximum  receiver  range 

KMN 

RAYSUM 

Index  of  minimum  range  of  all  NJ  ray  types 

KMNJ(J) 

RAYSUM 

Index  of  minimum  range  of  Jth  ray  type 

KMXJ(J) 

RAYSUM 

Index  of  maximum  range  of  Jth  ray  type 

NA 

INSERT 

No.  points  in  augmented  velocity  profile 

NB 

PLRAY 

No.  points  in  original  velocity  profile 

NCYC 

PLRAY 

Arrival  order  number 

NF 

PLRAY 

No.  frequencies 

NJ 

PLRAY 

No.  multipath  ray  types 

NL 

SECLIM 

No.  sectors  (less  one) 

NTH 

DIVSEC 

No.  rays  computed  in  a sector 

R(ITH,J) 

PLRAY 

Ray  range  (yd) 

RADE 

PLRAY 

Radius  of  earth  (ft) 

RCYC(ITH) 

INCR 

Ray  cycle  range  (yd) 

RMIN 

PLRAY 

Minimum  receiver  range  (yd) 

RP(ITH,J) 

PLRAY 

Range  derivative  (yd) 

RPCYC(ITH) 

INCR 

Derivative  of  RCYC(ITH)  (yd) 

RPSR(ITH) 

INCR 

Derivative  of  RSR(ITH)  (yd) 

RPUP(ITH) 

INCR 

Derivative  of  RUP(ITH)  (yd) 

RSR(ITH) 

INCR 

Range,  source  to  receiver  (yd) 

RTD 

PLRAY 

Conversion  factor,  radians  to  degrees 

RUP(ITH) 

INCR 

Range  of  portion  of  ray  cycle  above  receiver  depth  (yd) 

SIGN(J) 

PLRAY 

Algebraic  sign  of  ray  angle  at  true  source 

SLOSS (IF) 

SURF 

Surface  reflection  loss  (dB) 

SN 

INSERT 

Flag  indicating  whether  true  source  is  at  ZAA  or  ZBB 

TH(ITH) 

DIVSEC 

Ray  source  angle  (rad) 

TH1 

PLRAY. 

Ray  source  angle  at  inner  boundary  of  sector  (rad) 
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TH2 

PLRAY 

Ray  source  angle  at  outer  boundary  of  sector  (rad) 

VA(ITH) 

PLRAY 

Tangent  of  ray  source  angle 

VBM 

RAYSUM 

RCAUST 

SDUCT 

Tangent  of  source  angle  at  beam  location  (true  source) 

ZA(IA) 

INSERT 

Depth,  augmented  velocity  profile  (ft) 

ZAA 

INSERT 

Depth  of  ray  source  (ft) 

ZAD 

INSERT 

Bottom  of  surface  duct  (ft) 

ZB(IB) 

PLRAY 

Depth,  original  velocity  profile  (ft  or  m) 

ZBB 

INSERT 

Depth  of  ray  receiver  (ft) 

ZP  Cl B) 

CURFIT 

Parameter  of  V.P.  profile  segment 

ZR 

PLRAY 

True  receiver  depth  (ft) 

ZS 

PLRAY 

True  source  depth  (ft) 

COMMON  /DATA/ 

Symbol 

Definition 

DB(I , IBP , J) 

Bottom  loss  from  input  curves  (dB) 

DG (I , IBP , J) 

Bottom  grazing  angle  from  input  curves  (deg) 

FR(J) 

Bottom  loss  frequency  (Hz) 

COMMON  /COH/ 

Symbol 

Definition 

JCOH(IF) 

Coherence  control  parameter 

COMMON  /SDCT/ 

Symbol 

Definition 

EN 

No.  "trapped"  modes  in  surface  duct  (not  an  integer) 

FR 

Frequency  (Hz) 

THETA 

Phase  angle  in  depth  function  (rad) 

VO 

Difference  between  phase  velocity  of  perfectly 
trapped  first  mode  and  sound  speed  at  surface 

ZA 

Depth  of  bottom  of  surface  duct  (ft) 

INDEX  VARIABLES 

Symbol 

Application 

IA 

Augmented  velocity  profile 
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IB 

IF 

IL 

ITH 

J 

K 

KC 


Symbol 

CFAC 
FMAX 
HMIN 
I CASE 
I PLOT 
J1 

MAXK 

MJ 

NDCD 

NFR 

NREC 

RK 

RMAX 

S 

SGN 

T 

ZDMAX 

ZRO(I) 

ZSO 


Symbol 

CR 

CS 


Original  velocity  profile 
Frequency 

Sector  1 

Ray  in  sector 
Multipath  ray  type 
Range 

Range  relative  to  smallest  range  of  caustic  correction  interval 

PLRAY 

Definition 

Conversion  factor  (ft/m) 

Largest  specified  frequencies  (Hz) 

Intensity  corresponding  to  999  dB  propagation  loss 

Flag  indicating  whether  source  is  above  or  below  receiver 

Houston  plot  control  parameter 

Used  in  initiating  ray  ranges 

Maximum  permissible  no.  range  points  (251) 

Flag  (does  first  sector  extend  continuously  from  negative  to 
positive  source  angles?) 

End  code,  terminating  reading  of  V.P.  cards 
No.  frequencies  of  input  bottom  loss  curves 
No.  receiver  depths 
Receiver  range  (yd) 

Maximum  receiver  range  (yd) 

Salinity 

Algebraic  sign  of  ray  source  angle 
BT  temperature  (°F  or  °C) 

Max.  duct  depth  for  which  surface  duct  model  is  implemented  (1000  ft) 

True  receiver  depth  (ft) 

True  source  depth  (ft) 

INSERT 

Definition 

Sound  speed  at  true  receiver  depth  (ft/sec) 

Sound  speed  at  true  source  depth  (ft/sec) 
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Symbol 

G(IA) 

Symbol 

DEL 
DELO 
I MAX 
THT(I) 

Symbol 

CDT 

CUT 

DX(1A) 

DXP(IA) 

IDN 
I UP 

Q1.Q2.Q3 

S 

THBD 
W1,VV2 
VI,  V2 
ZDT 

ZUT 

Symbol 
AN  FAC 
ANPPC 
ANPPCA 
ANPPCB 
AO (J) ,A1(J) 
BLK 

BMF  (J  , JJ , I F) 


SECLIM 
Definition 
Sound  speed  gradient  (sec-1) 

DIVSEC 

Definition 

Ray  source  angle  increment  (rad) 

Dead  interval  at  sector  boundaries  (rad) 

Maximum  allowable  no.  rays  in  sector  (50) 

Temporary  array  of  source  angles  in  outer  portion  of  sector, 
symmetrically  located  with  respect  to  TH(I)  (rad) 

INCR 

Definition 

Temporary  storage  of  CA(IA)  at  bottom  of  layer  containing  lower 
ray  vertex  (ft/sec) 

Temporary  storage  of  CA(IA)  at  top  of  layer  containing  upper  ray 
vertex  (ft/sec) 

Increment  of  range  in  layer  IA  (yd) 

Derivative  of  DX  with  respect  to  VA  (yd) 

Index  of  layer  containing  lower  ray  vertex 
Index  of  layer  containing  upper  ray  vertex 
Constant  factors  in  formulas  for  DX  and  DXP 

Algebraic  sign  required  to  compute  ray  vertex  depth  from  vertex 
velocity 

Grazing  angle  at  bottom  (deg) 

Reciprocals  of  VI  and  V2  respectively 
Tangent  of  ray  angle  at  top/bottom  of  layer 

Temporary  storage  of  ZA(IA)  at  bottom  of  layer  containing  lower 
ray  vertex  (ft) 

Temporary  storage  of  ZA(IA)  at  top  of  layer  containing  upper  ray 
vertex  (ft) 

RAYSUM 

Definition 

Constant  factor  in  formula  for  ANPPC 

No.  range  points  per  cycle  of  ray  interference  pattern 

Value  of  ANPPC  at  ray  source  end 

Value  of  ANPPC  at  ray  receiver  end 

Parameters  in  intensity  computations  for  steep  angle  propagation 
Bottom  loss,  interpolated  to  Kth  range  point 
Beam  pressure  ratio  (equivalent  to  BMFAC(IF)) 
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COFACR 

COFACS 

CVK 

FAC 

FRACTA 

FRACTB 

FRACTR 

FRACTS 

HK 

HTRM(J , JJ , IF) 
ICOH 

ICTH(J) 

IND 
I SUM 
ITHMN 
ITHMX 
JA(J) 

JJ 

KMX 

KTH 

KTYP 

PHFAC 

PHIA 

PHIB 

PHIR 

PHIS 

RCOI1 

RK 

RMN 

RMNJ 

RMX 


Coherence  factor  at  receiver  end  of  ray 
Coherence  factor  at  source  end  of  ray 
Vertex  velocity,  interpolated  to  range  RK 

Attenuation  factor,  including  water,  surface,  bottom,  and  beam 
losses 

Absolute  value  of  FRACTS 

Coherence  attenuation  factor  computed  from  formula  involving  no. 
of  range  points  per  cycle 

Coherence  attenuation  factor  at  true  receiver  end  (applicable  when 
coherence  occurs  at  both  ends) 

Coherence  attenuation  factor  (including  effect  of  beam  pattern) 
Intensity  of  range  RK,  exclusive  of  coherence  factor 
Temporary  intensity  array 

Coherence  flag,  indicating  whether  coherence  occurs  at  neither  end 
(0),  ray  source  end  (1),  ray  receiver  end  (2),  or  both  ends  (3) 

Values  of  index  ITH  such  that  caustic  lies  between  ITH-1  and  ITH 

Parameter  indicating  which  of  multipath  rays  propagate  to  range  RK 

No.  of  multipath  rays  which  propagate  to  range  RK 

Minimum  value  of  ITH  in  ray  interpolation  loop 

Maximum  value  of  ITH  in  ray  interpolation  loop 

Identifier  indicating  which  of  the  multipath  ray  types  are  appli- 
cable in  coherence  computations 

Index  identifying  the  branches  of  range-VA  curve 

Largest  value  of  K among  all  NJ  multipath  ray  types  in  sector 

Same  as  ICTH(J) 

Flag  indicating  whether  special  formula  is  to  be  used  for  intensity 
computation 

Constant  factor  in  formula  for  phase  angle  in  interference  pattern 
Phase  angle  between  interfering  rays  at  ray  source  end  frad) 

Phase  angle  between  interfering  rays  at  ray  receiver  end  (rad) 

Phase  angle  between  interfering  rays  at  true  receiver  end  (rad) 

Phase  angle  between  interfering  rays  (rad)  • 

Nominal  limiting  value  of  horizontal  separation  of  interfering 
rays  at  source  or  receiver  depth  (ft) 

Receiver  range  at  index  K (yd) 

Minimum  range  among  all  NJ  multipath  ray  types  in  sector  (yd) 
Minimum  range  of  Jth  multipath  ray  type  in  sector  (yd) 

Maximum  range  among  all  NJ  multipath  ray  types  in  sector  (yd) 
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RMXJ 

RPK 

RPKK(J.JJ) 

SGN 

SLK 

THBD 

TRM 

VAK 

VAKK(J.JJ) 

VBK 

VO 

ZRCOH 

ZSCOH 


Symbol 
All 
CCL 
COE  FI 

COEF2 

CVC 

FAC 

HK 
I SGN 
K1 
K2 

K3 

K3T 

L 


1.2,  L3 

RC 

RCl 


Maximum  range  of  Jth  multipath  ray  type  in  sector  (yd) 

Interpolated  value  of  range  derivative  RP  at  range  RK 
Array  of  absolute  values  of  RPK 

Indicates  algebraic  sign  of  product  of  the  beam  functions  of  two 
interfering  rays 

Surface  loss  (dB) 

Grazing  angle  at  bottom  (deg) 

Intensity  in  final  form  for  insertion  into  H(K,IF)  array 
Interpolated  value  of  V A at  range  RK 
Array  of  values  of  VAK 

Tangent  of  angle  at  receiver  depth  of  interpolated  ray  to  RK 
Tangent  of  angle  at  surface  of  interpolated  ray  to  range  RK 
Horizontal  separation  of  two  interfering  rays  at  receiver  depth  (ft) 
Horizontal  separation  of  two  interfering  rays  at  source  depth  (ft) 

RCAUST 

Definition 

Square  of  Airy  function,  Ai(x) 

Constant  factor  2*(3tt)2^3 

Coefficient  relating  argument  of  Airy  function  to  horizontal 
distance  from  caustic 

Coefficient  relating  intensity  to  All 

Vertex  velocity  of  interpolated  ray  to  caustic  (ft/sec) 

Attenuation  factor,  including  water,  bottom,  surface,  and  beam 
losses 

Intensity 

Algebraic  sign  of  RPK 

Limiting  value  of  K on  shadow  side  of  caustic 

First  limiting  value  of  K on  insonified  side  of  caustic  (caustic 
correction  applies  unconditionally  between  K1  and  K2) 

Second  limiting  value  of  K on  insonified  side  of  caustic  (caustic 
correction  is  compared  with  ray  intensity  between  K2  and  K3) 

Cut-off  value  of  K3 

Range  index  beginning  with  1 at  K = K1  and  increasing  to 
|K3  - K1 | + 1 at  K = K3 

Values  of  L at  K = K2  and  K = K3 

Range  to  caustic  (yd) 

Limiting  range  on  shadow  side  of  caustic  (yd) 
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RC2 

RC3 

RK 

RPP 

RPPT 

TRM 

VAC 

VBC 

X 

XI 
X2 
X3 

Symbol 

ALPH1 

ALPH2 

ALPH3 

ALP1 

ALP  2 

CR 

CS 

FF 

FNIN 

HK 

HKK 

HM 

1ND 

RK 

TllliTR 
T1  lliTS 
U 
UK 


First  limiting  range  on  insonified  side  of  caustic  (yd) 

Second  limiting  range  on  insonified  side  of  caustic  (yd) 

Receiver  range  at  index  K (yd) 

Second  derivative  of  R with  respect  to  VA  (yd) 

Alternate  value  of  RPP  computed  from  adjacent  pair  of  rays 
Intensity 

Tangent  of  source  angle  of  interpolated  ray  to  caustic 

Tangent  of  receiver  angle  of  interpolated  ray  to  caustic 

Argument  of  Airy  Function 

Limiting  value  of  X on  shadow  side 

First  limiting  value  of  X on  insonified  side 

Second  limiting  value  of  X on  insonified  side 

SDUCT 

Definition 

Leakage  attenuation  coefficient  in  cylindrical  spreading  formula 
(dB/yd) 

Attenuation  coefficient  in  spherical  spreading  formula,  good  duct 
model  (dB/yd) 

Attenuation  coefficient  due  to  surface  roughness  (dB/yd) 

Overall  attenuation  coefficient,  cylindrical  spreading  formula 
(dB/yd) 

Overall  attenuation  coefficient,  spherical  spreading  formula, 
good  duct  model  (dB/yd) 

Constant  loss  term  associated  with  receiver  depth  (dB) 

Constant  loss  term  associated  with  source  depth  (dB) 

Frequency  (Hz)/100 

Minimum  frequency  for  which  surface  duct  computations  are  made  (Hz) 

Intensity  (exclusive  of  attenuation)  in  cylindrical  spreading 
formula 

Intensity  (exclusive  of  attenuation)  in  spherical  spreading  formula 
Minimum  intensity  at  which  surface  duct  computations  are  terminated 
Flag  to  avoid  repeated  computation  of  constants 
Range  at  index  K (yd) 

Phase  angle  in  receiver  depth  function  (rad) 

Phase  angle  in  source  depth  function  (rad) 

Product  of  Ui>  and  UR 
Receiver  depth  function 
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US 

ZA1 

Symbol 


Source  depth  function 

Depth  of  surface  duct  (ft/100) 

BEAM 

Definition 


BM(I , IF) 

NBF(IF) 
SINE  (I , IF) 
THETA(I) 


Beam  loss  factor  (pressure  ratio  relative  to  beam  axis)  for  the 
Ith  grazing  angle  and  IFth  frequency 

Number  of  points  required  to  specify  IFth  beam  pattern 

Sine  of  THETA(I)  for  IFth  frequency 

Grazing  angle  at  true  source  (deg) 
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VELOC.  Called  from  PLRAY.  Calculates  sound  speed  from  temperature  and 
salinity . 

CURFIT.  Called  from  PLRAY.  Fits  the  input  velocity  profile  points  with 
curvilinear  segments  which  join  with  continuous  slope.  PARAB,  FIT,  and  BRIDGE 
are  auxiliary  subroutines  called  by  CURFIT. 

INSERT.  Called  from  PLRAY.  Inserts  source  and  receiver  depths  into  velocity 
profile.  Of  these  two  depths  INSERT  selects  the  one  with  the  larger  sound  speed 
to  serve  as  the  source  for  ray  computations.  If  a surface  duct  is  present,  INSERT 
checks  the  source  and  receiver  depths  to  determine  whether  the  surface  duct  model 
is  to  be  implemented  and,  if  so,  sets  a flag. 

SECLIM.  Called  from  PLRAY.  Divides  the  angular  space  of  ray  source  angles 
into  sectors  bounded  by  limiting  rays  which  separate  rays  of  one  type  (e.g.,  RSR) 
from  rays  of  another  (e.g.,  bottom  bounce). 

D1VSEC . Called  from  PLRAY.  Computes  the  source  angles  of  a family  of 
non-uniformly  spaced  rays  in  each  sector. 

INCR.  Called  from  PLRAY.  Computes  range  and  range  derivative  increments 
for  each  ray  in  a sector.  The  increments  computed  are:  (1)  source  depth  to 
receiver  depth,  (2)  the  portion  of  a ray  cycle  above  receiver  depth,  and  (3)  a 
complete  ray  cycle. 

RAYSUM.  Called  from  PLRAY.  Processes  rays  of  a sector,  one  sector  at  a time. 
Searches  for  presence  of  caustics  and  calls  RCAUST  for  each  caustic  found. 

Computes  individual  ray  intensities  at  all  ranges  covered  by  rays  of  the  sector 
(exclusive  of  the  rays  covered  by  the  caustic  corrections) . Computes  the  in- 
coherent sum  of  multipath  ray  intensities  at  each  range.  When  applicable, 
computes  coherence  factors  and  resultant  semicoherent  intensities. 

RCAUST.  Called  from  RAYSUM.  Computes  intensities  over  a predetermined 
range  interval  on  either  side  of  a caustic. 

A12.  Called  from  RCAUST.  Computes  the  square  of  the  Airy  function  Ai(x), 
required  for  the  caustic  corrections. 

SURF.  Called  from  PLRAY.  Computes  the  surface  loss  as  a function  of  fre- 
quency and  sea  state.  (Does  not  include  any  dependence  on  ray  angle.) 

BOTTOM.  Called  from  INCR  and  RAYSUM.  Computes  bottom  reflection  loss  as 
a function  of  frequency,  ray  angle,  and  bottom  type.  Contains  bottom  loss 
curves  for  5 FNWC  bottom  types  internally. 

SOUCT.  Called  from  PLRAY.  Computes  intensities  from  a special  set  of 
formulas  in  lieu  of  conventional  ray  computations  when  either  source  or  re- 
ceiver depth  is  in  the  surface  duct  and  the  other  depth  does  not  exceed  1.8 
times  duct  depth. 

DEPTH.  Called  from  SDUCT.  Computes  depth  functions  required  for  surface 
duct  model. 
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BEAM. 

Part  1.  Called  from  PLRAY.  Reads  and  prints  beam  inputs. 

Part  2.  Called  from  RAYSUM,  RCAUST,  and  SDUCT.  Computes  beam 
pattern  functions. 

PLPLOT.  Called  from  PLRAY.  Plot  routine  for  Houston  plotter.  Plots 
propagation  loss  as  a function  of  range. 
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A complete  data  set  for  a normal  run  of  PLRAY  consists  of  8 cards  (or  card 
sets),  some  of  which  may  be  omitted,  depending  upon  the  options  selected  by  the 
user.  After  completing  the  computations  specified  by  the  first  data  set,  the 
program  returns  to  the  beginning  and  reads  the  first  card  of  a new  data  set.  If 
this  is  a blank  card,  the  program  stops.  However,  if  a complete  new  data  set 
is  supplied,  it  will  be  processed.  As  many  complete  data  sets  as  desired  may 
be  stacked  together  and  they  will  be  processed  sequentially.  A blank  card  must 
be  inserted  at  the  end  of  the  last  data  set  to  stop  the  program. 

The  instructions  which  follow  apply  to  a single  data  set.  The  first  input 
is  a card  containing  a set  of  control  integers.  The  second  input  is  a card  con- 
taining the  title  of  the  run  and  environmental  parameters  relating  to  the  ocean 
surface  and  bottom.  The  third  input  is  a set  of  cards  containing  externally 
supplied  bottom  loss  data.  The  fourth  input  is  a set  of  cards  containing  the 
sound  velocity  profile.  The  fifth  input  is  a card  containing  the  frequencies 
and  coherence  parameters.  The  sixth  input  is  one  or  two  cards  containing  the 
source  depth  and  receiver  depths.  The  seventh  input  is  a card  on  which  the 
desired  receiver  ranges  are  specified.  The  eighth  input  is  a set  of  cards  on 
which  the  source  beam  patterns  are  specified. 

Card  1,  control  integers 

NF,  NREC,  IPROF,  IUNIT,  IBEAM,  IPLOT,  IPR  Format  (1015) 

NF  = number  of  frequencies,  up  to  6 
NREC  = number  of  receiver  depths,  up  to  8 
IPROF  = environmental  input  control  parameter 

= -1,  VP  curve-fitting  only;  no  ray  computations 
= 0,  normal  operation 

= 1,  omit  both  VP  and  external  bottom  inputs;  bypass  Card  Sets  2-4; 

use  inputs  supplied  in  previous  data  set 
= 2,  omit  VP  inputs;  bypass  Card  Set  4;  use  VP  supplied  in  previous 

data  set 

= 3,  omit  external  bottom  inputs;  bypass  Card  Set  3;  use  external 

bottom  inputs  supplied  in  previous  data  set 
IUNIT  = unit  system  control  parameter 
= 0,  normal  operation 

= 1,  required  only  when  reading  temperatures  in  °C 
IBEAM  = beam  control  parameter 

= 0,  no  beam  patterns  involved;  bypass  Card  set  8 
= 1,  read  beam  patterns 

= 2,  bypass  Card  Set  8;  use  beam  patterns  supplied  in  previous  data  set 
IPLOT  = plot  control  parameter 
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=>  0,  no  plot 

= 1,  propagation  loss  plot  on  Houston  plotter 
IPR  = diagnostic  print  control 
- 0,  normal  operation 
= 1,  print  diagnostic  information 
= 2,  print  more  diagnostic  information 

Card  2,  environmental  header  card 

ISP,  IBP,  LAT,  IDR,  IDBT  Format  215,  10X,  A1 , 10X,  4A10) 

ISP  = surface  parameter  (sea  state) 

IBP  = bottom  parameter 

= -1,  infinite  bottom  loss  (100  dB) 

= 0,  zero  bottom  loss 
= 1 to  5,  FNWC  bottom  loss  curves 
= 6,  read  external  bottom  loss  curves 
LAT  = latitude  (deg),  a real  (decimal)  number.  Required  only  when  reading 
temperatures  and  salinities;  otherwise  only  for  identification.  De- 
fault value,  45  deg. 

IDR  = hemisphere  designation,  N or  S;  used  only  for  identification 
IDBT  = run  title,  up  to  40  alphanumeric  characters 

Card  Set  3,  external  bottom  inputs,  required  only  when  IBP  = 6 and  IPROF  = 0 or  2 

Card  3a,  frequencies  of  bottom  loss  curves;  must  be  specif iced  in  monotonically 
increasing  order 

NFR,  FR(J) , J = 1,  NFR)  Format  (15,  5X,  6F10.0) 

NFR  = number  of  frequencies,  up  to  6 
FR(J)  = frequency  (Hz) 

Card  Pair  3b,  bottom  loss  data 

(DB(I ,6,J)  I = 1,  15)  Format  (8F10.0) 

DB  (I,  6,  J)  = bottom  loss  (dB) . If  less  than  15  points  are  required  to 
specify  the  curve,  remaining  10-column  fields  may  be  left  blank;  2 
cards  required 

Card  Pair  3c 

(DG(I,6,J),  I * 1,  15)  Format  (8F10.0) 

DG(I,6,J)  = grazing  angle  (deg)  corresponding  to  loss  DB  (I , 6 , J ) . Last 
value  specified  must  be  90  degrees. 
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If  more  than  one  frequency  is  specified,  stack  cards  in  the  following  order: 

3b,  3c  for  first  frequency,  then  3b,  3c  for  second  frequency,  etc. 

Card  Set  4,  velocity  profile  deck,  one  card  for  each  point,  arranged  in  order  of 
increasing  depth. 

ZB(IB) , CB(IB) , T,  S,  NDCD  Format  (4F100.0,  II) 

ZB(IB)  = depth  (ft  or  m) 

CB(IB)  = sound  speed  (ft/sec  or  m/sec) 

T = temperature  (°F  or  °C) 

S = salinity 

NCDC  = end  code;  punch  a 1 in  column  41  of  last  card 
Inputs  may  be  either  sound  speeds  or  temperatures  and  salinities.  Use  same  set 
of  units  throughout,  either  metric  or  English.  When  specifying  temperatures  in 
°C,  I UNIT  must  be  1. 

Card  5,  frequencies  and  coherence  parameters. 

(F (IF) , IF  = 1,6),  (JCOH(IF) , IF  = 1,6)  Format  (6F10.0,  10X,  611) 

F ( I F)  = frequency  (Hz);  monotonic  order  not  required 
JCOH(IF)=  coherence  parameter 
- 0,  semicoherent 
= 1,  incoherent 

The  same  value  of  IF  applies  to  both  F(IF)  and  JCOH(IF).  Only  NF  values  of  each 
need  to  be  specified. 

Card  6,  source  and  receiver  depths 

(ZSO,  (ZRO(I) , I = 1,  NREC)  Format  (8F10.0 
ZSO  = source  depth  (ft) 

ZRO(I)  = receiver  depth  (ft).  If  8 receiver  depths  are  specified,  a second 
card  will  be  needed. 

Card  7,  ranges 

RMIN,  DR,  RMAX  Format  (8F10.0) 

RMIN  = minimum  range  (yd) 

DR  = range  increment  (yd) 

RMAX  = maximum  range  (yd) 

The  maximum  number  of  ranges  is  251.  If  more  than  251  are  specified,  RMAX  will 
be  automatically  adjusted. 

Card  Set  8,  beam  inputs 
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Card  8a 

(NBF(IF),  IF  = 1,  NF)  Format  (615) 

NBF(IF)  = number  of  points  required  to  specify  beam  for  IFth  frequency, 
up  to  100 

Card  Set  8b,  beam  ang les 

(THETA(I) , I = 1 ,N)  Format  (10F8.0) 

THETA(I)  = angle,  positive  upward.  The  angles  must  be  specified  in 
increasing  order  from  -90  to  +90  degrees 
N = NBF  (IF) 


Card  Set  8c,  beam  pattern  functions 

(BM(I,IF),  I = 1,  N)  Format  (10F8.0) 

BM(I,IF)  = pressure  ratio  relative  to  direction  of  maximum  response. 

When  more  than  one  beam  is  specified,  place  Card  Sets  8b  and  8c  for  the  first  beam 
immediately  after  Card  8a;  follow  these  with  Card  Sets  8b  and  8c  for  the  second 
beam,  then  the  third,  etc. 
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Let  zu  and  cu  refer  to  the  original  (uncorrected)  depth  and  sound  speed  in 
the  real  ocean,  zc  and  Cc  to  the  corrected  values  in  the  flat  ocean,  and  R to 
the  radius  of  the  earth.  As  described  in  reference  (a),  the  transformation 
consists  of  straightening  out  the  annular  ring  containing  the  ocean  into  a flat 
slab  in  such  a way  that  angles  are  locally  preserved.  The  ocean  surface  is 
unwrapped  without  changing  its  length.  Other  lines  of  constant  depth,  which  are 
arcs  of  concentric  circles,  must  be  stretched  while  being  straightened.  At  the 
same  time  they  must  be  spread  apart  in  order  to  preserve  local  directions.  The 
exact  correction  formulas  are 


z = -R  ln(l  - z /R) 
c u 


c = 


c 1 - z /R 
u 


(E-l) 

(E-2) 


Since  the  radius  of  the  earth  is  very  large  compared  with  the  ocean  depth, 
it  is  convenient  to  expand  the  logarithm  in  a series,  and  sufficiently  accurate 
to  throw  away  all  but  the  first  two  terms.  Thus, 

zc  = zu  (1  + Zu/2R)  (E-3) 

The  earth  curvature  correction  as  applied  in  PLRAY , consists  of  equations 
(E-3)  and  (E-2).  In  reference  (a),  equation  (E-2)  is  further  approximated,  but 
that  approximation  is  really  unnecessary. 
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STEEP  ANGLE  APPROXIMATION 
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The  fact  that  rays  propagate  in  approximately  straight-line  paths  at  steep 
angles  permits  the  use  of  an  approximation  which  avoids  the  necessity  of  comput- 
ing actual  ray  paths  and  thereby  saves  a significant  amount  of  computation.  At 
very  steep  angles,  where  the  ray  paths  are  very  close  to  straight  lines,  the 
horizontal  range  is  inversely  proportional  to  the  tangent  of  the  ran  angle.  At 
angles  beyond  about  30  degrees  or  so,  the  range  can  be  approximated  with  suffi- 
cient accuracy  by  a formula  of  the  form 


r = l/fa  v + a,) 
o a 1 


(F-l) 


where  aQ  and  aj  are  constants  and  va  is  the  tangent  of  the  ray  source  angle.  As 
applied  in  SUBROUTINE  RAYSUM,  r is  the  independently  specified  receiver  range, 
and  equation  (F-l)  is  solved  for  va  in  order  to  compute  the  range  derivative  r'. 
Thus, 


v = (1/r  - a )/a 
a 1 o 


(F-2) 


- — = - a r 
3v  o 

a 


(F-3) 


To  provide  smooth  transition  into  the  region  of  shallower  angles  where  exact 
ray  formulas  are  used,  the  constants  aQ  and  a^  are  evaluated  from  the  steepest 
ray  computed,  which  in  PLRAY  is  the  30-degree  ray.  Solving  equations  (F-l)  and 
(F-3)  for  aG  and  aj  yields 


ao  = ' r3o'/r30 


(F-4) 


ai  1/r30 


ao  V30 


(F-5) 


where  the  subscript  30  refers  to  the  30-degree  ray. 

The  above  approximation  is  one  of  the  excellent  novel  features  of  the  FACT 
model . 
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RAY  RANGES  AND  RANGE  DERIVATIVES 
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This  appendix  is  concerned  with  the  increments  in  range  Ar  and  range  deriva- 
tive Ar'  associated  with  the  passage  of  a ray  thourgh  a single  layer  of  the  veloc- 
ity profile.  In  PLRAY  the  independent  variable  with  respect  to  which  the  derivative 
is  taken  is  the  tangent  of  the  ray  source  angle,  denoted  by  va.  As  indicated  in 
the  text,  the  ray  source  may  be  either  the  true  source  or  the  true  receiver,  de- 
pending upon  which  has  the  larger  sound  speed. 

The  basic  formulas  for  the  range  and  range  derivative  are  derived  in  ref- 
erence (a),  pp.  32-39.  In  this  appendix  we  shall  apply  them  to  PLRAY  with  ap- 
propriate changes  in  notation.  The  relation  between  sound  speed  and  depth  in  a 
curvilinear  segment  (eq.  (64)  of  reference  (a)  will  be  written  in  the  form 


c -Jr^r r&T)z  <G-J> 

V p p 

where  c2,  gp,  and  zp  correspond  to  the  FORTRAN  symbols  CSP,  GP,  and  ZP.  In  a 
Type  I segment  all  tnree  of  these  parameters  are  positive;  in  a Type  II  segment 
gp  is  negative  and  Cp2  is  positive;  in  a Type  III  segment  gp  is  positive  and  Cp2 
is  negative.  (cp  itself,  which  would  be  imaginary,  has  no  physical  significance 
in  a Type  III  curve.) 

The  formula  for  Ar  in  a Type  I segment  (eq.  (73),  reference  (a)  is 

r = (sin'1  vl  - sin'1  w^/qj  (G-2) 


where  the  subscript  l refers  to  the  bottom  of  the  segment  and  u to  the  top,  and 


Equation  (G-2)  may  be  rewritten  in  a 


(G-3) 

(G-4a) 

(G-4b) 

form  more  suitable  for  computation 


r 


sin  1 (w«  /l-w  2 
£ u 


wu  ^ " wjiZ  )/<ll 


(G-S) 


The  formula  for  a Type  II  or  Type  III  segment  (eq.  (82)  of  reference  (a)  is 


r = 


(1/qj)  In 


qi(2&-zp)  + vj. 

q, (z  -z  ) + v 
M1  up  u 


(G-6) 


where  v and  v^  denote  the  tangents  of  the  ray  angle  at  the  top  and  bottom  of  the 
segment . 


The  formula  for  the  derivative,  applicable  to  all  segment  types  is  eq.  (75) 
of  reference  (a).  In  terms  of  the  ray  tangents  vy  and  v^,  it  is 
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where 


Z„  - Z Z - Z 

Au  “ Ar  ♦ q,  (4— E “ -V- E) 


2 2 2, 
q = c /(c  -c  ) 
3 p v p 


(G-7) 


(G-8) 


The  definition  of  Au,  eq.  (32)  of  reference  (a),  may  be  written  in  the  form 

d.Ar 


Au  = -cot  0 


a 30 


(G-9) 


where  @a  is  the  ray  source  angle.  The  derivative  is  taken  at  constant  depth. 
But 

. . 3Ar  3Ar  d0a  2„  3Ar 

= *'os  6a  00" 


3 v 39  dv 
a a a 


Hence 


And 


= - (c_2v  Au/c  2) 

3.  3^  V 

(G- 10) 

z.  -z  z -z  T 

q_  ( 4 E-  - -H_E  ) 

H3  v.  v ' 

SL  u -I 

(G-ll) 

At'  = -(ca2va/cv2)  [Ar  + **  ( 

In  SUBROUTINE  INCR  the  FORTRAN  symbols  for  Ar  and  Ar'  are  DX(IA)  and 
DXP(IA).  The  boundary  designated  with  the  subscript  l corresponds  to  the  FORTRAN 
index  IA+1,  and  the  subscript  u corresponds  to  IA.  The  subscript  a,  referring  to 
the  ray  source,  corresponds  to  the  symbol  A (or  AA)  in  the  FORTRAN. 
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APPENDIX  H 


intensity 


IN  THE  VICINITY  OF  A CAUSTIC 
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The  caustic  correction  formulas  contained  in  SUBROUTINE  RCAUSE  are  derived 
from  eq.  (83.43)  of  Brekhovskikh  (reference  (f)),  reproduced  below: 


2ll/6 


[tt/4  + 


w UQyj 


V(x) 


(H-I) 


where  ipj/  is  the  wave  function,  which  is  proportional  to  the  acoustic  pressure, 
w is  the  phase  of  the  wave,  obtained  from  the  WKB  approximation,  and  £ is  the 
horizontal  wave  number  which,  by  Snell's  law,  may  be  expressed  in  the  following 
forms : 


5 


0)  a <*>  r. 

— cos  0 = — COS  0, 

c a c,  b 

a b 


(H-2) 


In  eq.  (H-2)  u>  is  the  angular  frequency  (2irf)  and  the  other  symbols  have  the 
same  significance  as  in  appendix  I.  Other  symbols  in  eq.  (H-l)  are  defined  as 
follows : 


r = horizontal  range 


8 = — sin  0 = uv  /c 

O C a a V 

a 

(H-3) 

e = ^ sin  6b  = wvb/cv 

(H-4) 

W"(0  = 

(H-5) 

v(x)  = StT~  Ai(x) 

(H-6) 

where  Ai(x)  is  one  of  the  Airy  functions, 
is  given  by  eq.  (38.52)  of  Brekhovskikh 

The  argument  x of  the  Airy  function 

x-21'3 

which  may  be  written  in  the  form 

x = 21/3  (rc-r)/w'"U0)1/3 

(H-7) 

where  r^  is  the  range  to  the  caustic. 


Equation  (H-l)  is  based  on  a series  expansion  of  the  phase  w(£)  along  a 
horizontal  line  at  the  receiver  depth.  All  parameters  in  this  equation  except 
v(x)  are  evaluated  at  the  caustic.  In  particular,  £0  is  the  value  of  £ for  the 
ray  to  the  caustic.  The  argument  of  the  Airy  function  is  proportional  to  the 
horizontal  distance  from  the  caustic.  The  field  at  the  caustic  is  obtained  by 
inserting  x = 0.  Positive  values  of  x correspond  to  the  shadow  side  of  the 
caustic  and  negative  values  correspond  to  the  insonified  side. 

The  caustic  correction  intensity  is  unconditionally  computed  and  used  in 
lieu  of  the  uncorrected  ray  intensity  for  values  of  x between  3.5  and  -1.77. 

The  limit  3.5  lies  far  enough  down  on  the  exponential  decay  curve  of  the  Airy 


I 
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function  to  yield  intensities  of  negligible  value.  The  value  -1.77  corresponds 
to  a point  beyond  the  first  maximum  on  the  oscillatory  (insonified)  side  such 
that  under  ideal  conditions  the  intensity  would  be  equal  to  the  incoherent  sum 
of  the  intensities  of  the  two  rays  which  form  the  caustic,  thereby  assuring 
continuity  between  the  regions  of  the  ray  and  caustic  computations.  When  real- 
world  environments  are  inserted  in  PLRAY,  however,  the  conditions  for  the  series 
expansion  on  which  eq.  (H-l)  is  based  are  never  ideal,  and  a discontinuity 
usually  occurs.  To  avoid  the  discontinuity,  caustic  intensities  are  computed 
out  to  x = -2.2,  a point  well  down  toward  the  first  zero  of  the  Ai  function.  In 
the  interval  of  x from  -1.77  to  -2.2,  both  caustic  and  ray  intensities  are  com- 
puted and  at  each  range  the  larger  of  the  two  is  selected. 

To  convert  eq.  (H-l)  into  the  desired  intensity  for  PLRAY,  it  should  be 
noted  that  the  rms  acoustic  pressure  at  the  receiver  is 

pb  = i*n+i 

The  constant  factor  in  ip  + has  been  chosen  such  that  the  rms  acoustic  pressure  at 
the  reference  distance  or  1 yard  from  the  source  is 


In  normal  mode  solutions  the  intensity  ratio  is  considered  to  be  the  square  of 
the  ratio  of  these  two  pressures.  In  ray  solutions  a geometric  factor  involving 
cos  0 and  cos  0 enters  into  the  formulation,  leading  to  the  result 

3.  D 


H = 


Pb2/pCb  Capb2 
2.  " 2 
pa  /pca  cbpa 


(H-8) 


where  p is  the  density  of  water.  A discussion  of  the  ambiguity  of  these  two 
definitions  is  given  by  Pedersen  and  Gordon  in  reference  (g) . Since  the  current 
derivation  is  intended  for  use  in  a ray  program,  the  definition  given  in  eq.  (H-8) 
will  be  used. 


The  remainder  of  this  appendix  is  concerned  with  transforming  eq.  (H-l)  to 
the  form  in  which  it  is  programmed  in  SUBROUTINE  RCAUST.  At  the  outset  it  will 
be  noted  that  the  complex  exponential  factOT  may  be  dropped,  since  we  are  con- 
cerned only  with  its  magnitude.  Substitution  of  eqs.  (H-2)  to  (H-6)  into  eq . 
(H-l)  yields  , v / v 


= 2n'6A 


-1/3 


Ai(x) 


(H-9) 


It  may  be  shown  by  application  of  cq. 


(H-2),  together  with  the  definition 


that 


v = tan  0 
a a 


2 2 4 

“ Va  ca 


(H- 10) 


H-3 


J 


r f 
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where 


r»*  = 3 r 


(H-ll) 


Substitution  of  (H -10)  into  (H-9)  yields 


.♦  =2,2/v/6 


rv  v.c  ■ 
a b v 


The  intensity  ratio  , eq.  (H-8) , becomes 


. 4/3.1/3  3 

4rr  f ca 


rc,cv->v  v, 
b a b 


where 


VjL-Ca. 


2 

llLS a 
r" 


- 2/ 3 -2/ 3 / 2 

2tt  f ca  v,  c; 


Ai  (x) 


(i-  - rc) 


Ai(x) 


(H-12) 


(H-13) 


(H-14) 


In  eqs.  (H-13)  and  (H-14)  all  quantities  except  x and  r are  to  be  evaluated  for 
the  ray  which  propagates  to  the  caustic;  rc  is  the  range  to  the  caustic.  These 
equations  are  correct  in  their  present  form  if  all  lengths  are  in  yards.  How- 
ever, in  PLRAY  all  sound  speeds  are  in  ft/sec.  For  this  reason  the  right-hand 
side  of  eq.  (H-13)  must  be  multiplied  by  31/3  and  the  right-hand  side  of  eq. 
(H-14)  by  32/3. 


To  reexpress  eqs.  (H-13)  and  (H-14)  in  a form  similar  to  the  FORTRAN  coding 
in  RCAUST , let 

C = 2(3lQ2/3ca  / Va2ca\  1/3  (H-15) 

cv2  l r"  / 


cx  = f2/3C 


(H-16) 


F = CaCy 
1 3r  c.  v v. 
c b a b 

= C2FJf^/,3 


(H- 17) 


(H-18) 


H = C2Ai2(x) 


(H-18) 


x = Cj (r  - rc) 


(H- 19) 


The  similarity  of  the  symbols  should  facilitate  comparison  of  the  above  formulas 
with  the  FORTRAN  statements.  The  reader  is  reminded  that  the  subscript  a in  the 
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above  equation  corresponds  to  A or  AA  in  the  FORTRAN,  and  the  subscript  b 
corresponds  to  B or  BB. 

The  Airy  function  (Ai(x)  is  computed  from  the  standard  series  (reference  (h) ) 
for  values  of  x between  -2.2  and  1.0.  The  FUNCTION  AI2(X)  is  the  square  of  Ai(x). 
In  the  interval  from  1.0  to  3.5  a simple  approximation  has  been  developed  for  the 
computation  of  Ai^fx).  Let 


y = 6.709  + 3.04077  (x-2)  + 0.31288  (x-2)2  - 0.0195216  (x-2)3  (H-20) 


Then 


Ai2(x)  = e'y 


(H-21) 
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The  The  spreading  loss  factor  is  the  ratio  of  the  ray  intensity  at  the  receiver 
location  to  the  intensity  at  the  reference  point  one  yard  from  the  source,  ne- 
glecting all  losses  except  geometric  spreading.  PLRAY  uses  the  same  basic  formula 
as  the  old  NADC  Ray-Tracing  Program.  The  purpose  of  this  appendix  is  to  transform 
the  old  expression  contained  in  eq.  (38)  of  reference  (a)  into  the  form  appearing 
in  PLRAY.  The  spreading  loss  factor,  which  will  be  designated  by  the  symbol  H, 
is  the  reciprocal  of  the  argument  of  the  logarithm  in  eq.  (38).  With  appropriate 
changes  in  notation,  the  factor  is 


H = 


si  cos  6c 


r u sin  0 sin  0„ 
s R 


(1-1) 


where  9S  is  the  ray  angle  at  the  source,  0r  is  the  ray  angle  at  the  receiver, 
r is  the  horizontal  range,  si  is  the  length  of  one  yard,  expressed  in  the  same 
units  as  r and  u,  and  u is  defined  by  eq.  (32)  of  reference  (a). 


cos  fls  3r 
sin  0S  30s 


d-2) 


In  PLRAY  the  range  derivative  rs  taken  with  respect  to  v , the  tangent  of  0 

instead  of  0 itself  (0  in  reference  l) . a a 

a v s ' 


v = | tan  0 | 
a | a | 

Hence 


3r  _ 3r  dva  1 q 

30  3v  d0  cos2©-  av 
a a a a a 


(1-3) 

d-4) 


Let 


The  spreading  loss  factor  is  then 

2 3 . 

S]  cos  0a 

H r r'  sin  0^ 


2 3 n 

S1  cos  9a_ 
r r'  v,  cos  0, 
b b 


where 


Vb  = 


tan  Q,_ 


and  the  superscript  b refers  to  the  ray  receiver. 


d-5) 
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Finally,  applying  Snell's  law, 
ca  cb 


cos  0 


a 


cos  0,  cv 


where  ca  and  c^,  are  the  sound  speeds  at  the  ray  source  and  ray  receiver  depths 
and  c is  the  ray  vertex  velocity,  we  obtain 


H = 


si2  c 3 
S1  ca 


c z c,  v,  r r' 
v b b 


(1-6) 


Equation  (1-6)  is  in  the  form  in  which  it  is  programmed  in  SUBROUTINE  RAYSUM. 
Since  r and  r'  are  in  yards,  sj  has  the  value  unity.  The  FORTRAN  equivalents  of 
the  symbols  in  eq.  (1-6)  are: 


CAA  = c 


CBB  = c. 


CVK2  = c 2 
v 


VBK  = v. 


RK  = r 
RPK  = r' 
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COHERENCE  COMPUTATIONS 
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Like  the  FACT  model,  PLRAY  provides  the  option  of  semicoherent  ray  summation. 
The  significance  of  the  term  "semicoherent"  is  that  coherence  is  confined  to  the 
individual  rays  of  each  arrival  order;  the  resultant  intensities  of  the  different 
arrival  orders  are  added  incoherently.  The  semicoherent  approach  has  a certain 
amount  of  physical  reasonableness  in  that  the  various  rays  within  a given  arrival 
order  normally  propagate  by  paths  which  lie  close  together  and  are  very  different 
from  the  ray  paths  of  other  arrival  orders.  Another  good  feature  is  that  the 
semicoherent  ray  summation  preserves  the  smooth,  broad  features  of  the  interfer- 
ence pattern  which  frequently  can  be  correlated  with  experimental  data,  but  ignores 
the  rapid  fluctuations  which  in  the  real  world  have  only  a statistical  significance. 

The  most  accurate  way  to  compute  the  phase  difference  between  two  interfering 
rays  is  to  calculate  the  absolute  phases  from  the  ray  travel  times.  However,  both 
FACT  and  PLRAY  employ  a simple  approximate  method  which  does  not  require  the  cal- 
culation of  travel  time.  Coherence  is  computed  only  for  those  rays  which  are 
capable  of  striking  the  surface.  Consider  a pair  of  rays  whose  paths  to  the  re- 
ceiver differ  only  by  one  surface  reflection  at  the  source  end.  The  rays  leave 
the  source  at  approximately  equal  angles,  one  upward  and  the  other  downward.  The 
surface-reflected  ray  may  be  visualized  as  propagating  downward  from  an  image 
source  located  above  the  surface  of  the  water  symmetrically  with  respect  to  the 
real  source.  If  the  source  depth  is  a relatively  small  fraction  of  the  total 
depth  of  the  ocean,  the  two  ray  paths,  over  most  of  their  trajectories,  are  nearly 
identical  and,  in  the  immediate  vicinity  of  the  source  may  be  considered  as  ap- 
proximately parallel  straight  lines.  The  phase  difference  between  the  two  rays 
may  be  computed  to  a sufficient  approximation  from  the  difference  in  path  length, 
together  with  the  phase  reversal  resulting  from  the  surface  reflection.  Similar 
considerations  apply  to  the  case  in  which  the  surface  reflection  occurs  at  the 
receiver  end. 

If  the  source  (or  receiver)  depth  is  increased,  the  separation  between  the 
two  paths  becomes  greater,  and  the  coherence  between  the  two  rays  may  be  expected 
to  become  degraded.  In  FACT  and  PLRAY  the  effect  of  the  degradation  is  treated 
in  a somewhat  crude  manner.  The  horizontal  separation  of  the  two  rays  at  the 
source  (or  receiver)  depth  is  computed  and  compared  with  a pre-selected  threshold 
of  12,000  feet.  If  the  separation  does  not  exceed  the  threshold,  full  coherence 
is  assumed  to  exist.  Otherwise  the  rays  are  treated  incoherently. 

MULTIPATH  RAY  TYPE  DESIGNATIONS 

The  multipath  ray  types  are  identified  in  the  computer  program  by  the  index  J. 
In  the  zero  order  arrivals  there  are  two  rays  as  indicated  in  figure  Jl.  The 
direct  ray  is  identified  as  the  first  ray  (J  = 1)  and  the  surface-reflected  ray 
as  the  second  (J  = 2) . In  the  first  and  subsequent  order  arrivals  there  are  four 
rays.  As  illustrated  in  figure  J2  for  the  case  of  the  first  order  bottom-bounce 
arrivals,  the  first  ray  experiences  no  surface  reflections  at  either  end,  the 
second  reflects  at  the  source  end  (ZAA)  only,  the  third  reflects  at  the  receiver 
end  (ZBB)  only,  and  the  fourth  reflects  at  both  ends.  It  should  be  pointed  out 
that  the  special  case  (MJ  = 1),  in  which  there  is  only  one  ray  in  the  zero  order 
arrival  and  two  rays  in  the  first  and  subsequent  order  arrivals,  is  not  involved 
in  the  coherence  computations  because  in  that  case  the  rays  do  not  strike  the 
surface . 
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FIGURE  J2  - Multipath  Ray  Types,  First  Order  Arrivals. 
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Lven  though  the  total  number  of  rays  in  an  arrival  order  is  normally  two  for 
the  zeroth  order  and  four  for  the  subsequent  orders,  it  is  not  necessarily  true 
that  all  of  these  rays  are  actually  present  in  any  given  sector.  For  any  given 
source  depth  the  range  to  which  ray  1 propagates  is  always  less  than  the  ranges 
of  rays  2 and  3,  and  these  ranges  in  turn  are  always  less  than  that  of  ray  4. 

Hence,  at  the  inner  and  outer  edges  of  the  spread  of  ranges  covered  by  the  rays 
of  a sector,  there  is  usually  an  interval  where  one  or  more  of  the  rays  is  not 
present . 

It  is  also  not  necessarily  true  that  all  rays  present  in  a given  arrival 
order  are  actually  involved  in  the  coherence  computations.  The  reason  is  that 
one  or  more  of  the  rays  may  already  have  been  included  in  a caustic  correction. 

In  view  of  the  above,  it  is  necessary,  before  computing  the  coherence  factors, 
to  identify  the  individual  rays  involved.  Obviously  if  only  one  ray  is  present, 
there  is  no  coherence.  In  what  follows  we  shall  consider  the  cases  of  two,  three, 
and  four  rays  separately. 

TWO  RAYS 

Let  Ia  and  lb  be  the  individual  intensities  of  the  two  rays,  exclusive  of 
the  beam  pattern,  where  the  subscripts  a and  b refer  to  the  particular  pair  of 
rays  involved.  Let  ga  and  gb  be  the  respective  beam  pattern  functions  expressed 
as  pressure  ratios  relative  to  the  beam  axis.  It  is  assumed  that  the  symmetry 
of  the  beamis  such  that  ga  and  g^  are  positive  or  negative  real  numbers.  If  no 
beam  pattern  has  been  specified,  ga  and  g^  are  both  equal  to  unity. 

The  first  step  is  to  examine  the  horizontal  separation  of  the  rays  to  deter- 
mine whether  the  coherence  calculation  is  applicable.  The  horizontal  separation  is 

rsep  = 2z/tan  0o  (J-1} 

where  0O  is  the  ray  angle  at  the  surface  and  z is  the  depth  of  either  the  ray 
source  ZAA  or  the  ray  receiver  ZBB,  depending  upon  the  rays  involved.  ZAA 
applies  to  zero  order  arrivals  when  ZAA  is  less  than  ZBB  and  to  higher  order 
arrivals  when  the  rays  are  1 and  2 or  3 and  4.  ZBB  applies  to  zero  order  arrivals 
when  ZAA  is  greater  than  ZBB  and  to  higher  order  arrivals  when  the  rays  are  1 and  3 
or  2 and  4.  If  the  rays  are  1 and  4 or  2 and  3,  neither  pair  contains  a direct  and 
surface-reflected  ray  at  the  same  end,  and  it  is  assumed  that  no  coherence  is 
present.  Otherwise,  it  is  assumed  that  coherence  occurs  when  r is  less  than 
the  threshold  value.  seP 

If  coherence  occurs,  the  instantaneous  pressure  is 

P ■ /ra*a  - <Vb  <J-2' 

where  <f>  is  the  phase  difference  between  the  two  rays  due  to  the  path  difference 
and  the  minus  sign  accounts  for  the  phase  reversal  due  to  the  surface  reflection. 
Figure  J3  shows  the  geometry  of  the  parallel-line  approximation  from  which  <p  is 
computed.  S is  the  source  (or  receiver)  located  at  depth  z and  S'  is  its  image. 

The  direct  ray  is  SA  and  the  reflected  ray  is  SCD.  Since  the  actual  segment  SC 
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FIGURE  J3  - Approximate  Coherence  Model. 
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has  the  same  length  as  the  image  segment  S'C,  the  difference  in  path  length  S'B 
is  obtained  by  dropping  the  perpendicular  SB  on  the  image  ray.  The  phase  dif- 
ference is  found  to  be 

<P  = 4 irf  z sin  0 /c  (J-3) 

o o 

where  f is  the  frequency  and  cq  is  the  sound  speed  at  the  surface. 

The  coherent  intensity  is  given  by  the  square  of  the  complex  pressure, 
eq.  (J-2),  and  may  be  expressed  in  the  form 


= «a*a 


2 T 21 

* ‘b*b  > 


2 "Vb  Sagh 
I 2 2 + i p 2 
aKa  bKb 


cos  $ 


(J-4) 


Since  the  omnidirectional  intensities  of  the  component  rays  of  a given  arrival 
order  do  not  usually  differ  greatly,  we  shall  make  the  simplifying  assumption 
that  inside  the  parentheses  of  the  second  factor  of  eq.  (J-4)  they  may  be 
considered  equal.  With  this  assumption  the  resultant  intensity  becomes 


1 ' (Iaea2  * 'b*b2)  (‘  - 2 “s  *)  W-« 

' *a  * «b  ' 

The  first  factor  in  eq.  (J-5)  is  the  incoherent  sum  of  the  two  ray  intensities. 
The  effect  of  the  coherence  is  accounted  for  by  the  second  factor,  which  we  shall 
term  the  coherence  factor. 


THREE  RAYS 

In  PLRAY  the  three  ray  case  is  broken  down  into  a single  ray  and  a pair  of 
rays.  The  intensity  of  the  single  ray  is  added  incoherently  to  the  resultant 
intensity  of  the  pair.  In  all  possible  combinations  of  three  rays  it  is  found 
that  there  is  one  pair  of  rays  which  is  a candidate  for  coherence  at  the  ray 
source  end  and  another  pair  which  is  a candidate  for  coherence  at  the  ray  receiver 
end.  If  the  ray  source  depth  ZAA  is  less  than  the  ray  receiver  depth  ZBB,  the  co- 
herence is  calculated  at  the  source  end;  otherwise,  at  the  receiver  end.  The 
selection  rules  are  given  in  table  J-l.  The  ray  pair  is  then  treated  in  the  same 
manner  as  the  case  of  two  rays. 

TABLE  J-l 
RAY  SELECTION,  THREE-RAY  CASE 


Rays 

ZAA<ZBB 

ZAA>ZBB 

Present 

Pair 

Single 

Pair 

Singl 

1,2,3 

1,2 

3 

1,3 

2 

1,2,4 

1,2 

4 

2,4 

1 

1,3,4 

3,4 

1 

1,3 

4 

2,3,4 

3,4 

2 

2,4 

3 
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FOUR  RAYS 

Tests  art  made  of  the  applicability  of  coherence  at  both  the  ray  source  end 
and  the  ray  receiver  end.  If  the  coherence  is  not  to  be  computed  at  either  end, 
all  four  ray  intensities  are  summed  incoherently.  If  coherence  occurs  at  the 
ray  source  end  only,  rays  1 and  2 and  rays  3 and  4 are  treated  as  two  separate 
pairs.  If  coherence  occurs  at  the  receiver  end  only,  rays  1 and  3 and  rays  2 
and  4 are  treated  as  separate  pairs. 

In  the  general  four  ray  case  the  instantaneous  pressure  is 


P = *1  gi  ‘ g2  el<^A  ' *3  g3  el<f>B  + *4  g4 


(J-6) 


where  and  $3  are  given  by  eq.  (J-3)  with  z equal  to  ZAA  and  ZBB  respectively. 
The  resultant  intensity,  which  is  the  square  of  the  modulus  of  p,  may  be  written 
in  the  form 


I = I, 


I- 


2 1 1 2 glg2  + *77  g3g4) 


cos  *. 


2(^77glg3  + *77g2g4} 

17 ' 

2(A2i3  g2g3  /iLi4  gtg4) 

ri 

2^  g2g3  + *ir4  glg4 


COS  <f 


B 


COS.  if  COS  <f 


B 


sin  <t>A  sin  4>R 


(J-7) 


where  1^  is  the  incoherent  sum 

II  = Ilgl2  + I2g22  + I3g32  + I4g42 


(J-8) 


In  order  to  reduce  eq.  (J-7)  to  a more  tractible  form,  it  is  necessary  to 
make  some  simplifying  approximations.  In  PLRAY  the  beam  pattern  is  applied  to 
the  true  source.  The  true  source  may  be  either  the  ray  source  or  the  ray  re- 
ceiver, depending  upon  which  the  two  depths  has  the  larger  sound  speed.  If 
the  beam  is  at  the  ray  source,  then  rays  1 and  3 have  negative  source  angles, 
while  rays  2 and  4 have  positive  source  angles.  Let 


g+  = h ( g2  * Sj 


(J-9a) 


&ah  (gx  + g3) 


(J-9b) 


Since  all  four  rays  have  very  nearly  the  same  trajectory,  it  is  reasonable  to  make 
the  approximations 


g2 .S  g4  S g+ 

gl  £*3  *«- 


(J-lOa) 

(J-lOb) 
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With  these  approximations  eq.  (J-7)  becomes 

1 = !i  (2 " cos  *a)  (1  ■ cos  v (j-n) 

If  the  beam  is  at  the  ray  receiver,  let 

g+  = *s  (g3  + g4)  (J-12a) 

g.  - *i  (gj  ♦ g2)  (J-  12b) 

If  we  make  the  assumptions 

g3  % g4  % g+  (J-13a) 

% *2  % 8-  (J- 13b) 

Equation  (J-7)  simplifies  to 


1 - li  ^ - g2f2?g  2 cos  ^b)  (1  " COS  V (J_14) 

EFFECT  OF  INADEQUATE  SAMPLING 

If  the  specified  receiver  range  points  at  which  the  propagation  loss  is  to 
be  computed  are  too  far  apart  relative  to  the  cycles  of  oscillation  of  the  ray 
interference  pattern,  there  exists  a sampling  problem.  If  the  sampling  is  to  be 
adequate,  there  should  be  at  least  a half  dozen  or  so  range  points  per  cycle.* 

The  problem  of  inadequate  sampling  is  handled  arbitrarily  by  cutting  down  the 
amplitude  of  the  interference  oscillations  by  multiplying  each  of  the  cosine 
functions  by  an  attenuation  factor  F which  varies  from  0 to  1.  In  the  FACT 
model  F is  zero  when  the  number  of  points  per  cycle  is  less  than  8/3;  it  is  1 
when  the  number  of  points  per  cycle  is  greater  than  6 and  varies  linearly  for 
intermediate  values.  In  PLRAY  the  piecewise  linear  function  is  replaced  with  a 
continuous  inverse  hyperbolic  function 

_ T,  1.2(4.333-n  )1  ,T  ir, 

F = 1/  1 + e ppc'  (J -15) 

where  nppC  is  the  number  of  receiver  range  points  per  cycle.  A plot  of  this 
function  is  shown  in  figure  J4. 

It  now  remains  to  develop  a suitable  formula  for  nppC.  The  number  of  cycles 
is  obtained  by  dividing  the  phase  angle  $ in  eq.  (J-3)  by  2 it.  Thus, 

nc=  2fz  sin  <|>o/co  (J-16) 


♦Recent  experience  has  suggested  that  this  number  may  be  higher  than  necessary. 

J-9 
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FIGURE  J4  - Attenuation  Factor  as  a Function  of  the  Number  of  Range  Points 
Per  Cycle  of  Ray  Interference  Pattern. 
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where  z is  either  the  ray  source  depth  ZAA  or  the  ray  receiver  depth  ZBB,  as 
applicable.  The  number  of  cycles  per  unit  range  is 


dn 

dn 

dv 

1 dn 

c 

c 

a 

c 

dr 

dv 

dr 

r ' dv 

a a 


(J- 17) 


where  va  is  the  tangent  of  the  ray  angle  at  the  ray  source,  and  r'  is  the 
derivative  of  the  range  with  respect  to  v„.  From  eq.  (J-16)  and  Snell's  law 
it  may  be  shown  that 

2 2 

dn  2fzsine  cos&  0 2fzv  c 
c _ a a o _ a a 

dv  c sinQ  v c 3 l.J-18) 

a o o o v 

where  0a  and  ca  are  the  ray  angle  and  sound  speed  at  the  ray  source,  vQ  is  the 
tangent  of  0O,  and  cv  is  the  ray  vertex  velocity.  If  the  range  increment  be- 
tween successive  receiver  locations  is  Ar,  the  number  of  points  per  cycle  is 
found  to  be 

The  attenuation  factor  F is  applied  directly  as  a multiplier  of  cos  <|>g 
eq.  (J-ll)  and  cos  <£«  in  eq.  (J-14) . When  applied  to  the  terms  involving  the 
coefficient  2gagj,/ (g^+gfc2)  in  eq.  (J-5)  and  the  coefficient  2g+g_/ (g+2+g_2) 
in  eqs.  (J-ll)  and  J-14),  the  attenuation  factor  F is  tested  against  the  absolute 
value  of  the  coefficient,  and  the  smaller  of  the  two  is  selected.  If  the  factor 
F is  chosen,  it  is  given  the  algebraic  sign  of  the  coefficient. 


J-ll 
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The  surface  duct  model  incorporated  in  PLRAY  is  divided  into  two  parts,  a 
poor  duct  model  and  a good  duct  model.  The  poor  duct  model  is  applicable  to 
those  combinations  of  frequency  and  duct  depth  for  which  the  propagation  is 
dominated  (except  at  very  short  ranges)  by  a single  normal  mode.  The  good  duct 
model  applies  to  cases  where  significant  contributions  at  all  ranges  of  interest 
are  made  by  several  normal  modes.  The  cirterion  for  determining  when  to  use  one 
model  or  the  other  is  based  on  a parameter  which  we  have  termed  the  number  of 
trapped  modes.  Before  describing  the  two  models,  it  will  be  advantageous  to 
develop  the  concept  and  derive  the  formula  for  the  number  of  trapped  modes. 

NUMBER  OF  TRAPPED  MODES 


The  surface  duct  model  is  derived  from  the  same  two- layer  bottomless  ocean 
as  forms  the  basis  of  the  Pedersen-Gordon  normal  mode  program  (reference  (i)) 
and  of  a similar  locally  generated  NAVAIRDEVCEN  normal  mode  program.  The  first 
layer  is  the  surface  duct,  whose  thickness  is  za.  The  sound  speed  and  gradient 
at  the  surface  are  cQ  and  . The  second  speed  varies  with  depth  according  to  the 


.Z)  \ c.  / 


(K-l) 


The  sound  speed  ca  at  the  bottom  of  the  duct  is  found  by  substituting  za  into 
eq.  (K-l).  The  sound  speed  in  the  second  layer  varies  in  a similar  manner,  the 
gradient  y^  being  negative. 

The  mode  depth  functions  are  solutions  of  the  depth  equation 


d2u  + / ai 
dz2  \ c(z) 


2 - k2y  u = 0 


(K-2) 


where  k is  a parameter  which  takes  on  a set  of  discrete  values  (eigenvalues) 
corresponding  to  the  normal  modes.  With  the  substitution 


ZU)  - ij  (k2  - ) 

a \ c„  / 


(K-3) 


where 


, 2 ,1/3. 

q = (2u>  Yx)  /cQ 


(K-4) 


equation  (K-2)  transforms  to  Stokes'  equation 


- Zu  = 0 


whose  solution  is  a linear  combination  of  the  Airy  function  Ai(Z)  and  Bi(Z), 


u(Z)  = A Ai (Z)  + B Bi (Z) 


K-2 
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If  it  is  assumed  (as  an  approximation)  that  the  duct  is  thick  enough  that  u(Z) 
decays  to  a negligible  value  at  the  bottom  (Z  = Za) , then  the  coefficient  B 
vanishes  and  u(Z)  consists  only  of  the  Ai  function. 

To  determine  the  discrete  modes,  we  apply  the  boundary  condition  at  the 
surface,  u(0)  = 0,  (equivalent  to  neglecting  the  impedance  of  the  air),  or 

Ai  [Z(0)]  = 0 (K-5) 

The  exact  solution  for  the  first  mode,  which  will  of  interest  later,  is 

Z(0)  = -2.3381074  (K-6) 

At  this  point  we  shall  be  concerned  with  an  approximate  analytic  solution 
obtained  from  the  asymptotic  approximation 


where 


Ai(Z)  -v  - - 1 — si"  (C  + */4) 

A (-z)1/4 


3/2 


(K-7) 

(K-8) 


The  boundary  condition,  eq.  (K-5),  is  satisfied  if  the 
eq.  (K-7)  has  the  value  nir , where  n is  an  integer- -the 
stitution  of  the  formula  (K-3)  for  Z,  the  condition  of 

t0  be  3 / \ 

k2)  ■ (n  - 1/4)* 


argument  of  the  sine  in 
mode  number.  Upon  sub- 
mode  resonance  is  found 


(K-9) 


Expressing  k in  terms  of  the  mode  phase  velocity  cv  (the  vertex  velocity  of  the 
equivalent  ray)  permits  the  further  transformation 


(n  - H)v 


(K- 10) 


Equations  (K-9)  and  (K-10)  define  discrete  values  of  k and  cv,  corresponding 
to  the  discrete  modes  of  a normal  mode  solution.  For  the  present  purpose,  how- 
ever, a useful  measure  of  the  quality  of  the  duct  can  be  obtained  by  setting  cv 
equal  to  the  sound  speed  at  the  bottom  of  the  duct  and  solving  eq.  (K-10)  for 
a nonintegral  value  of  n.  The  result  is 


n 


,5/2  1/2 


( K- 11) 


Inserting  representative  values  yj  = 0.016  sec'1,  c0  = 5000  ft/sec,  and  express- 
ing the  frequency  and  duct  depth  in  terms  of  the  reduced  valuables 


f'  = f/100  and  z ' = z /100 
a a 


(K-12) 


3/2 


wc  obtain 


n 


0.25  ♦ 0.067462  f'z  » 
a 


(K-13) 
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The  parameter  n may  be  loosely  termed  the  number  of  "trapped"  modes. 

Actually,  of  course,  none  of  the  modes  are  truly  trapped.  There  is  a continuous 
transition  region  between  the  lower  order  modes  (if  any)  which  are  tightly 
trapped  and  the  higher  order  modes  which  are  very  leaky.  The  value  of  n given 
above  lies  somewhere  near  the  middle  of  the  transition  and  as  such  is  a reasonable 
measure  of  the  number  of  modes  which  contribute  significantly  to  the  resultant 
intensity. 

Runs  made  on  the  NAVAIRDEVCEN  surface  duct  normal  mode  program  have  indicated 
that  a value  of  1.5  for  the  number  of  trapped  modes  constitutes  a reasonable  divid 
ing  line  between  the  poor  duct  region  and  the  good  duct  region. 

POOR  DUCT  MODEL 


The  poor  duct  model  is  a semi-empirical  simplification  of  the  normal  mode 
solution  for  the  first  mode.  The  formula  for  the  pressure  can  be  recognized  as 
having  the  form  of  the  normal  mode  solution.  Several  of  the  key  parameters  have 
been  evaluated  as  functions  of  frequency  and  duct  depth  by  empirical  curve-fitting 
to  the  results  of  runs  on  the  NAVAIRDEVCEN  normal  mode  program. 


The  mode  wave  function,  which  is  proportional  to  the  acoustic  pressure,  has 
the  form  


* =>/wu(zS)  u(zR]  fil 


(ir/4  - kr) 


(K- 14) 


where  u(z  ) and  u(z^)  are  functions  of  the  source  and  receiver  depths  z and  z . 
To  account  for  the  mode  attenuation  due  to  leakage  of  energy  out  of  the  auct,  the 
wave  number  k is  made  a complex  number 


k = w/c  - ia 
v l 


(K- 15) 


For  propagation  within  the  surface  duct,  the  phase  velocity  c may  be  replaced  by 
Cq  without  undue  loss  of  accuracy. 

The  mode  intensity  is  equal  to  the  square  of  the  modulus  of  the  wave  function 
In  forming  the  modulus  it  should  be  noted  that 

I.  I 2.2  2 . 

Ikl  ^ w / c0  +«1  ~ w/c0 


i(ir/4  - kr) 


-a.  r 
-el 


It  should  also  be  noted  that  the  mode  is  attenuated  by  absorption  of  energy  in  the 
water  «w'r  and  scattering  at  the  surface  as ' r in  addition  to  the  diffraction  leak- 
age ajr  mentioned  above.  Inserting  these  extra  terms  and  converting  the  attenua- 
tion coefficients  from  nepers/yd  to  db/yd,  we  obtain  for  the  intensity 


where 


«,  - ci 

otj  = mode  leakage  attenuation  coefficient, 

K-4 


(K- 16) 
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ag  = surface  scattering  attenuation  coefficient, 

a = water  attenuation  coefficient 
w 

all  in  db/yd. 

The  water  attenuation  coefficient  is  needed  elsewhere  in  the  program  and  is 

in  COMMON.  The  surface  scattering  attenuation  coefficient  is  taken  from  the  FACT 

model.  Allowing  for  the  difference  in  the  units  of  r between  FACT  and  PLRAY,  the 

formula  for  a is 
s 


W /f/1000  z. 


C K - 17) 


where 


0.00444,  sea  state  3 
W = ( 0.00666,  sea  state  = 3 
0.00888,  sea  state  3 J 


(K- 18) 


The  remainder  of  this  section  is  devoted  to  evaluation  of  the  depth  function 
u(z)  and  the  leakage  attentuation  coefficient  a . The  depth  function  u(z)  is  a 
modification  of  the  WKB  approximate  formula 


where 

and 

Now, 


l 


U(2)  -T-J 


const . 

fz~zn 

v ' %7 


1/4 


sin 


ft 

3 Yi 
0)  1 


c = phase  velocity  of  first  mode 


c(z) 


_L  _L_  - Cclf  + c)(cv  - c) 
2 " 2 2 2 


3/2 


(K- 19) 


For  propagation  in  the  surface  duct. 


c 'v  c 'v  c_ 

v 'v  'X/  0 


and 


_£XL 


/I  - 2y . 


'V. 


Therefore 


K-5 
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2 3^ 


J_  % 2(cv-c0) 

2 3 

c c_ 

v 0 


With  these  approximations  eq.  (K-19)  becomes 


u(z)  = - 


const . 


(WY1Z)% 


(£]  Si"|fc)5/2^[(VC°,3/2  ' (CVC»-T>2)W2Ji 

' ( ' ' (K-20) 


For  a further  transformation,  let 


c - c 
v 0 


(K-21) 


The  turning  point  of  the  mode  depth  function  occurs  at  the  depth  z T such  that 


Ac  - Y1zt  = 0 
or 

zT  = Ac/yj 

With  these  substitutions  the  depth  function  becomes 

, , const.  , / c0  V . ( 25/2tt 

u(z)  = ' (T-z/'z,  sln  < V 3/2 


(K-22) 


sil 


2S/2, 


f*  V2  , . ,3/2 

fAc  l-(i-z/zT) 


CK-23) 


Inserting  the  numerical  values  yi  = 0.016  sec  and  Cg  = 5000  ft/sec  and  the 
reduced  variables  f'  and  z ' (eq.  (K-12)),  we  obtain 

SL 


u(z)  = 


(l-z/z_ 


sin  < 0.10472  f' Ac“ 


-(l-z/zT)' 


(K-24) 


where  the  coefficient  A includes  all  the  constant  factors. 


Equation  (K-24)  was  used  as  a starting  point  for  the  development  of  a suitable 
function  for  the  poor  duct  model.  To  serve  as  a guide  in  the  development,  a set  of 
runs  were  made  on  the  NAVAIRDEVCEN  surface  duct  normal  mode  program  for  a variety 
of  frequencies  and  duct  depths.  The  previously  stated  values  were  assumed  for  Cq 
and  yj.  The  value  assumed  for  Y2  was  -0.1  sec-*.  Runs  were  made  for  a set  of 
source  depths  within  the  duct  and  for  a set  of  receiver  depths  extending  to  1.8  z . 
Computations  were  made  not  only  of  propagation  loss,  but  also  of  the  mode  depth 
functions. 

The  parameters  A and  Ac  were  evaluated  as  functions  of  frequency  and  duct 
depth  by  comparison  with  the  normal  mode  depth  functions,  using  suitable  curve- 
fitting techniques.  In  addition,  eq.  (K-24)  was  modified  in  the  vicinity  of  the 
turning  point  (where  the  WKB  approximation  diverges  to  infinity)  and  at  depths 
below  the  turning  point. 

K-6 
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Evaluation  of  Ac.  The  parameter  Ac  is  expressed  as  the  sum  of  two  terms 


where 


Ac  = ACq  + Ac^ 


(K-2S) 


Ac 


0 


c 


V oo 


C 

o 


(K-26) 


AC1  = cv  " cv»  (*-27) 

and  cVoo  is  the  phase  velocity  which  would  exist  if  the  mode  were  completely  trapped 
(no  leakage).  We  may  evaluate  Acq  from  theoretical  considerations.  From  eq.  (K-S) 
the  condition  of  resonance  for  the  first  mode  is 


c 2 ln 

f.2/3 

i i \ 

where 

0 \y 

:> 

“T 2 1 

C0  cv»  J 

Zq  = -Z(0)  = 2.3381074 
Equation  (K-28)  solves  to 


which,  upon  insertion  of  the  appropriate  values,  becomes 


(K-28) 


(K-29) 


Ac 


0 


5000 


I 1 - 

l A - 0.0692/f^-5 


(K-30) 


The  parameter  Acj  was  evaluated  by  an  empirical  fit  to  the  normal  mode 
eigenvalues. 


where 


a 


11 


Acx  = eUll  - ai2n) 

2.9  + 0.885  z • - 0.087  z ’2 
a a 

4.328  + 0.18  z ' , 
a 


z ' < 4 
a 


z >4 
a — 


(K-31) 


a 


12 


5.5 


and  n is  given  by  eq.  (K-13). 


Evaluation  of  coefficient  A. 


A = a2i  ♦ e’(a22  + a23 


K-7 
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where 


. = t 

0.04  + 0.39  f'  - 0.006  f'2. 

f' 

< 

2.5 

21  1 

1 

0.075  ♦ 0.01  f’. 

f' 

> 

2.5 

a = 

[ 1.415  - 0. 31f ' , 

f' 

< 

2.0 

22 

0.795 

f' 

> 

2.0 

f 0.1808  + 0.3996  f ’ , 

f' 

< 

2.0 

a25 = ; 

' 0.594  + 0.193  f'  , 

f' 

> 

2.0 

Formulas  for  the  Depth  Function.  Three  different  modifications  of  eq.  (K-24) 
are  used,  depending  upon  the  value  of  the  depth  z.  Let 


0.8  z 


Zj  = smaller  of  •{ 


z2  = 1.2  za 


a 


0.9  z„ 


Region  1,  z <_  z^ 


u(z*)  = ^ sin  <fr(z) 


where 


f Cl  - z/zT) 


h 


0 82 


D = larger  of  j 

<J>(z)  = 0.10472  f'Ac3/2  £l  - (1  - z/zt)3/2J 
Region  2,  z ^ < z <_ 


Let 

Then 

where 


Uj  = ufZj)' 

u(z)  = Uj  ea31^Z  Z\)^i 


a31  = 0.244  - 2.0  n + 4.842  e 


-3.58  n 


Let 


Region  3,  z_  < z <_  1.8  z 

* 3 


u2  = u(z2) 


(K-33) 

(K-34) 

(K-35) 


(K-36) 

(K-37) 

(K-38) 


(K-39) 
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u( z)  = u2  ea32(z  ' Z2)/Za 


(K-40) 


where 


a^2  = -0.444  + 3.42  e 


-3.49  n 


Evaluation  of  a].  The  result  of  the  curve-fitting  for  the  leakage  attenua- 


tion coefficient  is 


where 


ctj  = e(a41  + a42)/1000 


. , . ,0.1586 

a.  = 3.7  / z ' 

41  a 


(K-41) 


-3.237  n + 0.7965  e 


2.78  n - 10.422  n 


Alternate  Formula  for  Poor  Duct.  The  definition  of  a poor  duct  in  terms  of 
the  contribution  of  a single  mode  tends  to  break  down  at  short  ranges,  where 
higher  order  modes  become  important.  These  modes  have  little  or  no  effect  at  the 
longer  ranges  because  they  are  highly  attenuated.  From  a ray  point  of  view  we  are 
concerned  with  the  region  of  spherical  spreading  and  the  transition  to  cylindrical 
spreading.  A reasonable  approximation  to  the  intensity  in  this  region  is  provided 
by  the  following  formula 


H2  = —~2  sin2<J>(zs)  sin2<J>(zR)  • 10'0-1^al  + as+  °w-*r 
r 


(K-42) 


where  <J>(z)  is  the  phase  angle  of  eq.  (K-36) . Equation  (K-42)  is  applied  at  ranges 
less  than  2500  yd.  At  ranges  beyond  2500  yd,  it  is  multiplied  by  the  factor 


1Q-(r  - 2500)/5000 


(K-43) 


Resultant  Intensity.  At  each  receiver  range  the  two  intensities  Hj  and  H2 
are  compared  and  the  larger  of  the  two  is  selected.  Because  of  the  tendency 
toward  high  mode  attenuation  (01)  in  the  cylindrical  spreading  formula  (K-16), 
the  attenuation  factor  (K-43)  is  included  in  the  spherical  spreading  formula 
(K-42)  to  prevent  H2  from  becoming  dominant  at  the  longer  ranges.  When  the  range 
reaches  a value  such  that  the  intensity  is  less  than  10-^^,  no  further  computa- 
tions are  made  in  the  surface  duct. 

GOOD  DUCT  MODEL 

Since  no  simple  theoretical  formulation  appears  to  be  feasible  for  the  case 
where  many  modes  are  involved,  the  good  duct  model  was  formulated  by  empirical 
curve-fitting  to  a pair  of  formulas,  one  based  on  cylindrical  spreading  and  the 
other  on  spherical  spreading.  The  two  formulas  are 


1 
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„ 1 Tn-O.lfCo.  + a + a )r  + Kc  + K„1  .2,  .2, 

H,  = — • 10  1 1 s w S RJsin  <J>  sin  4>„ 

1 r s K 


u 2 -0.1(a_  + a + a )r  . 2 , .2 

H2  = — j • 10  2 s w'  sin  <J>s  sin  <J>R 

r 


(K-44) 


(K-45) 


where  a. , a , and  a are  the  same  as  in  the  poor  duct  model  and 
is  , w 


= 0.10472  f’  ac03/2  fl  - (1  - V2T0)3/^] 


4>r  = 0.10472  £•  Acq3/2  |~L  - (1  - VtT0)3/2J 


subject  to  the  limitation  that  neither  of  these  angles  may  exceed  tt/2  radians,  and 
^0  = A C0^Y1 

To  obtain  formulas  for  the  parameters  cu,  K„,  and  K , we  define  the  auxiliary 
relations  R 


:(z)  = e'6-0  (Z/Za  ’ °-92> 


xs  = x(zs) 

XR  = x(zR} 

^ = smaller  of  x_  and  x 

u K 


1.5 

a2  " 1000  (1  + x^) 

K(z)  = 52-/-  13  x(z) 
v 1 + x(z) 


K(zs)  ■ 


Kr  -Key 

As  in  the  case  of  the  poor  duct.  Hi  and  H2  are  compared  at  each  receiver  range 
and  the  larger  of  the  two  is  selected.  Computations  in  the  duct  are  terminated 
when  the  intensity  drops  below  the  minimum  value  of  10~13. 

CUT-OFF  FREQUENCY 

If  the  surface  duct  is  so  leaky  that  it  fails  to  trap  any  significant  amount 
of  energy,  there  is  no  point  in  bothering  with  the  computations.  In  the  PLRAY 


A 
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surface  duct  model  the  criterion  for  bypassing  duct  computations  is  expressed  in 
terms  of  a minimum  cut-off  frequency  fmin»  below  which  the  surface  duct  model  is 
bypassed.  The  cut-off  frequency  is  expressed  as  a function  of  the  reduced  duct 
depth  za'  (eq.  (K-12))  by  the  following  formula 

fmin  = 200/2 a'2'4  (Hz)  (*-46) 

The  above  relationship  corresponds  roughly  to  a mode  attenuation  of  about  IS 
dB/kyd.  Figure  K1  is  a plot  of  eq.  (K-46) . 


K-ll 


FIGURE  K1  - Cut-off  Frequency  as  a Function  of  Duct  Depth. 
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AP2  is  a normal  mode  program  which  accepts  an  arbitrarily  specified  ve- 
locity profile  in  the  same  form  as  PLRAY  and  FACT.  The  bottom  may  be  speci- 
fied either  in  the  form  of  a set  of  up  to  ten  homogeneous  fluid  type  layers  or 
in  the  form  of  the  same  type  of  empirical  bottom  loss  curves  as  are  used  in 
PLRAY  and  FACT. 

The  basic  approach  to  the  normal  mode  solution  is  essentially  the  same  as 
that  of  the  NAVAIRDEVCEN  Two-Layer  and  Three-Layer  Programs,  reference  (j). 

The  integration  in  the  complex  plane  is  based  on  the  same  branch  cut,  with  the 
result  that  the  solution  is  obtained  in  the  form  of  the  sum  of  an  infinite 
number  of  discrete  modes  and  a branch  line  integral  whose  value  under  most 
practical  operating  conditions  is  negligible.  As  in  reference  (j),  the  input 
velocity  profile  is  fitted  with  segments  in  which  the  square  of  the  index  of 
refraction  varies  linearly  with  depth,  so  that  the  solution  of  the  depth  equa- 
tion is  obtained  in  terms  of  Airy  functions.  However,  instead  of  attempting  a 
complete  algebraic  formulation  of  the  type  employed  in  reference  (j),  which 
obviously  becomes  impractical  with  more  than  three  la>ers,  the  solution  is  ob- 
tained by  an  iterative  procedure  similar  to  that  of  the  old  Program  for  Arbi- 
trary Velocity  Profiles,  reference  (k) . The  approach  bears  considerable  simi- 
larity to  that  of  the  discrete  mode  portion  of  D.  C.  Stickler’s  normal  mode 
program,  reference  (1).  The  multi-layered  and  empirical  bottoms  are  treated 
in  essentially  the  same  manner  as  in  reference  (j). 

AP2  is  dimensioned  for  500  modes.  To  avoid  unnecessary  computation  and 
to  make  space  available  for  higher  order  modes  above  the  500th,  AP2  contains 
an  automatic  procedure  for  bypassing  those  lower  order  modes  whose  excitation 
is  too  weak  to  contribute  significantly  to  the  resultant  intensity.  The  pro- 
cedure applies  to  those  modes  whose  phase  velocities  are  sufficiently  smaller 
than  the  sound  speeds  at  the  source  depth  and  all  specified  receiver  depths  to 
produce  the  required  amplitude  attenuation  beyond  the  mode  turning  point.  The 
mode  computations  are  automatically  terminated  at  the  first  mode  whose  phase 
velocity  exceeds  a value  equal  to  10  times  the  sound  speed  at  the  bottom  of 
the  water  or  after  the  500th  mode  has  been  computed,  whichever  occurs  first. 

When  checked  against  the  programs  of  references  (j)  and  (k) , AP2  has  been 
found  to  give  identical  outputs  for  the  same  inputs. 

Within  the  500-mode  limitation  AP2  has  the  same  versatility  as  ray  pro- 
grams such  as  PLRAY  and  FACT.  It  is  completely  automatic  and  does  not  contain 
any  special  inputs  which  require  the  user  to  have  a knowledge  of  normal  mode 
theory  or  of  program  idiosyncrasies.  It  is  planned  to  document  AP2  with  a 
formal  report  at  the  earliest  opportunity. 
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